Rudolf Zeller

Linear-scaling total-energy calculations with the tight-binding Korringa-Kohn-Rostoker Green function method

Keywords: electronic structure, metallic systems user supplied): complex energies, linear scaling, iterative solution electronic structure, linear scaling, iterative solution, complex energies

A complex-energy broadening scheme with quartic dependence on the broadening parameter is presented for Green-function density-functional electronic-structure calculations. The scheme is applied in a recently developed linear-scaling algorithm based on the tight-binding Korringa-Kohn-Rostoker Green function method where it leads to considerable computational saving compared to the standard Fermi-Dirac broadening with quadratic dependence. The linearscaling algorithm is applied for total-energy and spin-moment calculations for large supercells and the advantages of quartic broadening scheme are discussed.

Introduction

Standard density-functional calculations are usually restricted to about a thousand atoms even on the best supercomputers available today, since the computational effort increases with the third power of the number of atoms in the system. In the last years a number of codes have appeared [1][2][3][4] which overcome this bottleneck and achieve linear-scaling effort by utilizing the decay of the density matrix [5], which is exponential in semiconducting and insulating materials. For metallic systems, where (at zero temperature) the decay is only algebraical, a linear-scaling algorithm based on the tight-binding (TB) Korringa-Kohn-Rostoker (KKR) Green function method has been suggested recently [6]. This algorithm exploits the sparsity of the TB-KKR matrix, uses an iterative solution technique for the KKR matrix equations and applies a spatial truncation of the Green function to achieve linear scaling. For large Cu and Pd supercells as model systems it has been found [6] that the quasi-minimal residual (QMR) method [7,8] leads to convergent iterations, if the charge density is obtained by complex-energy integration [9] with non-zero electronic temperature [10], and that the total-energy error introduced by the spatial truncation of the Green function is in the meV range, if the truncation region contains about a few thousand atoms.

Here the previous study, which was limited to the minimal choice l max = 2 for the angular-momentum cut-off, is extended to l max = 3, which is more reasonable for total-energy calculations, and to a spin-polarized material to illustrate that the algorithm is also useful for the calculation of spin moments. Since the main computational effort for the iterative solution appears at the complex energy with smallest imaginary part, which is proportional to temperature T , it is important that T can be chosen as large as acceptable for the desired total-energy accuracy. Therefore, it is one aim of this paper to present and explain a complex-energy integration scheme, where the usual T 2 behaviour of the calculated energy is replaced by a T 4 behaviour, and to show that this scheme considerably reduces the computational cost for the iterative solution.

Basic equations

In the full-potential KKR Green function method space is divided into nonoverlapping space-filling cells around the atomic positions R n . In each cell the Green function can be written [11] in a multiple-scattering representation as

G(r + R n , r + R n ; E) = δ nn G n s (r, r ; E) + LL R n L (r; E)G nn LL (E)R n L (r ; E), (1)
where r and r are cell centred coordinates and L stands for the angular-momentum indices l and m and r and r are cell centred coordinates. G n s (r, r ; E) and R n L (r; E) denote single-scattering Green function and wavefunctions which depend only on the potential V (r) inside cell n and can be calculated with linear-scaling effort by solving single-site integral equations as described in [11,12] or in [START_REF] Zabloudil | Electron Scattering in Solid Matter: A Theoretical and Computational Treatise[END_REF]. For spinpolarized potentials the Green function and wavefunctions carry an additional spin index which, however, will be suppressed in the following to simplify the notation. The density n(r + R n ) inside cell n is obtained [10,11] by integration

n(r + R n ) = - 1 π Im ∞ -∞ dEf (E -E F , T )G(r + R n , r + R n ; E + i ), (2)
where

f (E -E F , T ) = (1 + exp((E -E F )/kT )) -1
is the Fermi-Dirac function for temperature T and E F denotes the Fermi energy. The infinitesimal quantity → 0+ guarantees that the singularities of the Green function on the real energy axis are avoided and that the Green function in the correct energy half plane is used. The Green function matrix elements G nn LL (E) are obtained by solving the matrix equation

G nn LL (E) = G r,nn LL (E) + n L L G r,nn LL (E)∆t n L L (E)G n n L L (E), (3)
where G r,nn LL (E) are the Green function matrix elements of suitably chosen reference system and ∆t n LL (E) are t-matrix differences which can be calculated from the potential and the reference potential, both restricted to cell n. The system (3) of linear equations represents an algebraic Dyson equation with complex, non-Hermitian matrices of dimension N (l max + 1) 2 .

The standard reference system in the KKR method is free space. Here the Green function matrix elements G 0,nn LL (E), traditionally called structure constants, are known analytically, but decay unfavourably slowly with distance between site n and n . Because of the slow decay all elements in the KKR matrix t -1 -G 0 have to be taken into account and because of that iterative solution of (3) is much more expensive than direct solution since both require a computational effort proportional to N 3 , which must be repeated in each iteration, whereas it is needed only once for direct solution. A useful reference system [14] with rapidly decaying structure constants in the energy range relevant for density-functional calculations consists of an infinite array of repulsive potentials (here chosen with a constant height of 8 Ry), which are confined to nonoverlapping muffin-tin spheres around the sites R n . The Green function matrix elements of the reference system, also called screened structure constants, can be calculated in real space with linear scaling effort by solving the Dyson equation

G r,nn LL (E) = G 0,nn LL (E) + n L L G 0,nn LL (E)t r,n L L (E)G r,n n L L (E) (4)
separately for each site n by use of a finite cluster of N cl repulsive potentials with the t matrices t r,n L L (E). For the present study the choice N cl = 13 was used which includes the central site and its nearest neighbours. Due to the rapid decay of G r,n n L L (E) with distance |R n -R n | other sites n were neglected in (4). This neglect reduces the dimension of the linear system (4) to N cl (l max + 1) 2 , leads to a sparse TB-KKR matrix with sparsity degree N cl /N and makes the computational effort proportional to N it N cl N 2 for the iterative solution of (3) compared to N 3 for the direct solution so that essentially N 2 scaling is achieved by iterative solution as has been discussed already by Smirnov and Johnson [15].

An important feature of iterative solutions is that for each site n the Green function matrix elements can be determined independently so that iterative solutions are ideally suited for massively parallel computing. A straightforward iteration of (3) corresponds to the Born iteration of scattering theory which often diverges, in particular for energies with small imaginary part. This problem can be avoided if one of the known iterative solution techniques for systems of linear equations is used which works for non-Hermitian matrices. One such technique, which was used successfully in the current investigation, is the quasi-minimal residual (QMR) method [7,8] in its transpose free version, which was used also by Smirnov and Johnson [15].

In order to arrive at a linear-scaling algorithm the principle of nearsightedness of electronic matter [START_REF] Prodan | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Kohn | [END_REF], which is exploited in other linear-scaling methods [1][2][3][4][5] for the density matrix, was used for the Green function G(r, r ; E) which decays at the complex energy mesh points in a similar fashion as the finite-temperature single-particle density matrix ρ(r, r , T ) as can be seen from the equation

ρ(r, r , T ) = - 1 π Im ∞ -∞ dEf (E -E F , T )G(r, r ; E). (5)
If this integral is evaluated by using the integrand at the Matsubara energies and on the complex-energy contour as described in the appendix, the decay of ρ(r, r , T ) is mainly determined by the decay of G(r, r ; E F +iπkT ), where E F +iπkT is the first Matsubara energy with smallest imaginary part, since the Green function decays faster on all other integration mesh points E i which have larger values of ImE i . Thus a neglect of the Green function in (2) for large distances |rr | corresponds to a neglect of the finite-temperature density matrix for similar distances. Since the single-scattering wavefunctions in (1) are only multiplicative factors, a truncation of the Green function directly corresponds to a neglect of Green function matrix elements G nn LL beyond a chosen distance d cut , which means that in (3) only O(N tr N ) elements G nn LL appear instead of O(N 2 ) so that the computational effort is reduced by a factor N tr /N . Here N tr is the number of atoms which are included in the truncation region given by |R n -R n | < d cut . The total effort for all N atoms is then proportional to N it N cl N tr N . This increases linearly with N since N it approaches a constant value for large N as will be illustrated in the next section. It should be pointed out that because of the singularities of the Green function on the real energy axis the use of complex energies is unavoidable and that the required number N it of iterations increases with decreasing imaginary part of E i .

Total-energy results and iteration behaviour

Previously [6] the proposed linear-scaling algorithm has been studied for large Cu and Pd supercells for the minimal choice l max = 2 of the angular-momentum cut-off. Here similarly large supercells were studied with l max = 3, which is more relevant for total-energy calculations. This represents an enormous reduction of the computational effort in the present model study compared to realistic systems with inequivalent atoms. According to the concept of equivalent k point meshes [18] the self-consistent potentials for the supercells agree with the ones of the simple cubic unit cells with four atoms if 5984 k points in the irreducible part of the simple cubic Brillouin zone are used. This means that self-consistent potentials for the large supercells could be obtained in an rather inexpensive way. An important issue for the proposed algorithm is how fast the iterations converge. The main computational work consists in matrix-vector multiplications, which were repeated for each site n and angular-momentum component L independently until the solution is obtained with the prescribed precision, here specified by the residual norm ||r|| = 10 -6 for the QMR iterations. The number N it of iterations depends on system size and on the imaginary part of the energy E i . The highest number of iterations is required at the first Matsubara energy E F + iπkT . In the previous investigation [6] for l max = 2 it has been found that the dependence of N it on the number N tr of atoms in the truncation region could be fitted to an exponential function of the form

N it (N tr ) = N ∞ it -α exp(γN 1/3 tr ) (6)
with three temperature dependent parameters N ∞ it , α and γ. In the present study for l max = 3 the same behaviour was found for T = 800 K and T = 1600 K. For truncation regions with more than 20000 atoms, the results of N it (N tr ) for T = 400 K scattered too much so that the fits became ill-defined and reliable fit parameters could not be extracted. Note that the values for N ∞ it given in Table 1, which correspond to averages over the (l max + 1) 2 = 16 angular-momentum components, are larger than the values obtained previously for l max = 2. However, those values, which stood for the number of QMR iterations, cannot be compared directly to the ones in Table 1, since here the more meaningful numbers of matrixvector multiplications are displayed which are about twice as large as the numbers of calls to the QMR subroutines displayed previously. An important issue for the proposed linear-scaling algorithm concerns the question how accurate the total energy can be calculated if spatial truncation for the Green function is applied. Figure 1 shows the total-energy error as function of the number N tr of atoms in the truncation region and illustrates that the error is smaller than 2 meV if truncation regions with a few thousand atoms are used and that the proposed algorithm thus can be expected to be useful for linearscaling total-energy calculations for large metallic systems. Another effect seen in figure 1 is that, for truncation regions with up to 10000 atoms, higher electronic temperature does not lead to substantially lower total-energy errors. Probably the zero-temperature algebraical decay of the Green function (and density matrix) dominates the additional exponential decay caused by temperature up to truncation distances of approximately 10 times the face-centred-cubic lattice constant.

http
It should be mentioned that it was not necessary to make self-consistent calculation for the total-energy error caused by the truncation because the modified energy functional Ẽ[n(r)] = E[n(r)]-E F n(r) -N el was used instead of the standard energy functional E[n(r)]. This functional is extremal for density variations, even if the total number N el of electrons changes, and has been demonstrated to yield accurate results in impurity calculations [19], where charge-neutrality deviations of about 0.01 to 0.1 electrons exist.

Quartic broadening scheme

From the results in Table 1 it is clear that the computational effort, which is proportional to the number of iterations (matrix-vector multiplications), is considerably reduced for higher electronic temperature. It has been argued [20] that large temperature can be detrimental for total-energy calculations, in particular for magnetic systems, where magnetism vanishes at higher temperature. Therefore, if one is interested in the total energy E tot = E(0), it is advantageous to replace the broadening function f (E -E F , T ) in (2), which leads to the temperature expansion

E(T ) = E tot + α 2 T 2 + α 4 T 4 + ... (7)
(see e. g. [START_REF] Ashcroft | Solid State Physics: Advances in Research and Applications[END_REF]) by a more powerful broadening scheme without T 2 term. This has been recognized in the past and various schemes have been suggested [START_REF] Methfessel | [END_REF][23][24]. Another, perhaps obvious idea, which was investigated here, is to calculate the energy E(T ) at two different temperatures and to eliminate the T 2 term in ( 7) by subtraction, for instance as

4 3 E(T ) - 1 3 E(2T ) = E tot -4α 4 T 4 -... . (8)
It is of course impractical to subtract the two large numbers on the left hand side of ( 8) and the calculation for the right hand side of ( 8) must be implemented directly in the charge density integration (2) by replacing f (E -E F , T ) with 4f (E -E F , T )/3 -f (E -E F , 2T )/3 in the integrand. On the real energy axis this difference introduces negative occupancies as the Methfessel-Paxton scheme [START_REF] Methfessel | [END_REF], which are avoided in the cold-smearing T 3 scheme of Marzari et al. [23]. For the complex-energy integration over the Green function the negative values occurring in the difference of two Fermi-Dirac functions are practically unimportant as will be explained in the appendix. In this way a useful T 4 scheme is obtained which also avoids cumbersome entropy calculation needed if E tot is approximated [24] by adding energy and free energy and dividing by two. The advantage of ( 7) over ( 8) is illustrated in figure 2, where the differences ∆ 2 = E tot -E(T ) and ∆ 4 = 4E(T )/3 -E(2T )/3 -E tot for several temperatures are plotted for face-centred-cubic Cu and Ni. The lines in figure 2 are fits to the functions E tot + α 2 T 2 or E tot + α 4 T 4 with two fit parameters each, E tot and α 2 or E tot and α 4 . Here it is important to note that enough k points and enough mesh points on the complex energy contours were used. Otherwise the calculated total energies are not accurate enough to obtain reasonable fits, which extrapolate the T 2 and T 4 behaviour to almost the same values for E(0) = E tot .

Figure 2 clearly shows that (8) leads to considerably better results than (7), for instance, the values for T = 1600 K obtained by (8) are nearer to E tot than the values for T = 800 obtained by (7). This is of importance for the computational effort Solid and open squares are for 800 K calculated with (7) and 1600 K calculated with (8). The lines, which connect the results for T = 800 K, serve as guide for the eye.

in the proposed linear-scaling algorithm since the necessary number of iterations is considerably smaller if T can be doubled (see Table 1). An interesting issue concerns the question whether the advantage of ( 8) over ( 7) appears also for other physical quantities, for instance for the spin moments in magnetic materials. This question was studied here for face-centred-cubic Ni. Whereas total-energy calculation did not require self-consistent results because of the modified functional Ẽ[n(r)] used, self-consistency was necessary for calculation of spin moments. The spin moments calculated with broadening functions f (E -E F , T ) and 4f (E -E F , T )/3 -f (E -E F , 2T )/3 are plotted in figure 3a. This figure shows that spin moments moments are reasonably well described for not too large temperature, in particular if the T 4 scheme is used. (Note that the Curie temperature T c of Ni is not determined by vanishing spin moments, but by disorder of the moments so that temperatures larger than T c can be used). Figure 3b shows that the additional error introduced by the truncation of the Green function required for the linear-scaling algorithm is smaller than the temperature error so that the linear-scaling algorithm can be expected to be useful also for magnetic properties of large systems if the truncation region contains enough atoms.

Discussion and summary

It was shown for the example of large supercells that iterative solution of the KKR Dyson equation and a spatial truncation of the Green function lead to an algorithm useful for linear-scaling total-energy calculations for large metallic systems. The requirements for the algorithm are the sparsity of TB-KKR matrix, the principle of nearsightedness of electronic matter utilized by the truncation, iterative solution, and the use of finite-temperature complex-energy integration to obtain the charge density. The behaviour of the algorithm was illustrated for model systems consisting of large Cu, Pd and Ni supercells which could be particularly easily studied because all atoms in the supercells were equivalent and because an appropriate choice of k points in the Brillouin zone facilitated the determination of the self-consistent potential for the supercells. It was shown that total-energy and spinmoment errors introduced by truncation can be smaller than 2 meV and 10 the exponential dependence on E, both functions have high amplitude oscillations for complex energies, which, however, do not exist on the contour for the Fermi-Dirac function because of its periodicity f (E + 2JiπkT, T ) = f (E, T ) so that the Fermi-Dirac function on the contour is real and as well behaving as on the real energy axis.

For the numerical integration the contour interval from E 0 + 2JiπkT to infinity was divided into several subintervals. On part of these subintervals the Fermi-Dirac function could be replaced by 1 and standard Gauss-Legendre integration rules could be used, whereas for other subintervals special Gauss rules had to be constructed for the weight functions w 1 (E) = f (E -E F , T ) and w 2 (E) = 4f (E -E F , T )/3 -f (E F -E, 2T )/3. For w 1 (E) one subinterval was chosen from E F + 2JiπkT to infinity and for w 2 (E) from E F + x 0 + 2JiπkT to infinity. Here x 0 = 2 ln(2 + √ 7)kT ≈ 3.071906kT is the only point on the contour, where w 2 (E) changes sign, so that the construction of Gauss rules was not complicated by sign changes of the weight function. Thus the correct choice of subintervals in the complex-energy plane with x 0 being at the boundary of two subintervals completely removes the problem of negative values of the broadening function which leads to the trouble of negative occupancies in basis-set methods [23]. The construction of Gauss integration rules is straightforward (see e. g. [START_REF] Gautschi | Orthogonal Polynomials: Computation and Approximation[END_REF]), although numerical steps might be susceptible to rounding errors because of illconditioning. Here the rules were constructed according to the following prescription. Monic orthogonal polynomials p n (x) = x n + ... satisfy a three term recurrence p n+1 (x) = (x -α n )p n (x) -β n p n-1 (x) starting with p -1 = 0 and p 0 = 1, where the coefficients are given by α

n = C -1 n xp 2 n (x)w(x)dx and β n = C -1 n-1 C n with C n defined as C n = p 2
n (x)w(x)dx and β 0 = w(x)dx. For an N point integration rule the mesh points x i are given by the roots of the polynomial p N (x) and the integration weights are determined by

w i = C n-1 p N (x i )p N -1 (x i ) , (A1)
where the prime denotes derivative with respect to x. Numerically the values of x i and w i were obtained by a computer algebra system with 90 digits precision to avoid rounding error problems. In total the number of Gauss points on the contour was about 50, which kept the total-energy errors arising from energy integration below 1 µRy.

Figure 1 .

 1 Figure 1. Total-energy error per atom as function of the number of atoms contained in the truncation region. Solid and open squares are for T = 800 K and 1600 K, diamonds for T = 400 K. The lines, which connect the results for T = 800 K, serve as guide for the eye.

Figure 2 .

 2 Figure 2. Energy difference ∆ 2 = Etot -E(T ) and ∆ 4 = 4/3E(T ) -1/3E(2T ) -Etot per atom as function of temperature. Solid and open squares are for ∆ 2 and ∆ 4 and the lines are fitted to a T 2 and T 4 behaviour.

Figure 3 .

 3 Figure 3. (a) Moment difference M (T ) -M (0) as function of temperature. Solid and open squares are for the T 2 nd T 4 schemes. The values up to 1200 K were used to obtain the lines by fitting a T 2 and T 4 behaviour. (b) Moment error per atom as function of the number of atoms contained in the truncation region.Solid and open squares are for 800 K calculated with(7) and 1600 K calculated with(8). The lines, which connect the results for T = 800 K, serve as guide for the eye.

Table 1 .

 1 Parameters determined by fitting the number of iterations required to obtain a relative residual norm ||r|| < 10 -6 at the Matsubara energy EF + iπkT to a function of the form(6) for the temperatures T = 800 K and T = 1600 K.

		N ∞ it (800)	α(800)	γ(800)	N ∞ it (1600)	α(1600)	γ(1600)
	Ni -	7568	7679	-0.0066	694	781	-0.0506
	Ni +	1151	1292	-0.0401	518	597	-0.0591
	Cu	1347	1630	-0.0537	519	707	-0.0985
	Pd	2339	2563	-0.0271	566	673	-0.0661

  The supercells contained 4 × 323 = 131072 atoms arranged in the appropriate face-centred-cubic geometry and were obtained by repeating a simple cubic unit with four atoms 32 times in the three space directions. A single point k = (1/4, 1/4, 1/4) × 2π/a was used in the irreducible part of the Brillouin zone. The supercell lattice constant a was chosen as 11.276, 11.568 and 124.74 nm for Ni, Cu and Pd, which is 32 times the experimental lattice constant of these elements. Since all atoms in the supercells are equivalent, the iterative solution of the Dyson equation (3) is needed for only one value of n .

Page 7 of 15 http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

  -3 µ B if the truncation region contains about a few thousand atoms. It would be desirable to extend the present model study to realistic systems with inequivalent atoms, where obviously much more computing power is needed, which is expected to be available on massively parallel supercomputers.

	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

http://mc.

manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

Appendix A. Construction of complex-energy integration rules

Instead of integrating (2) along the real energy axis, the integral can be evaluated in the complex energy plane because the Green function is an analytical function except for singularities for real energies. For the present study the contour starts on the negative real energy axis at energy E 0 below the valence and above the core states and thus includes only contributions of valence states. Contributions of core states were calculated separately by treating these states in an atomic like fashion. From E 0 the contour goes parallel to the imaginary axis up to E 0 + 2JiπkT , where J was a chosen small even integer. On this part of the contour the Fermi-Dirac function was replaced by the constant value 1 since for the temperatures used the quantity exp((E 0 -E F )/(2kT )) was smaller than 10 -20 . From E 0 + 2JiπkT the contour goes parallel to the real axis to infinity. Compared to the integral from E 0 to infinity along the real axis, the integral on the contour differs by the contribution of the residues arising from the poles of the Fermi-Dirac function at the first J Matsubara energies E j = E F +(2j-1)iπkT with j = 1, 2, ..., J, which are added separately. Compared to other broadening functions [START_REF] Methfessel | [END_REF]23] the Fermi-Dirac function is more advantageous in the complex-energy plane despite the necessity to include pole contributions, which do not exist if a Gauss function is used. Due to


Uli Buchenau
email: buchenau-juelich@t-online.de

Herbert R Schober

An atomic mechanism for the boson peak in metallic glasses

Keywords: boson peak, glass transition, metallic glasses, vibrational properties boson peak, glass transition, metallic glasses PACS numbers: 64.70.Pf, 77.22.Gm

published or not. The documents may come

An atomic mechanism for the boson peak in metallic glasses

As yet there is no generally accepted explanation of the boson peak in the neutron or Raman scattering intensities of glasses. The boson peak is a broad peak at an energy transfer of a few meV, where simple crystals show only sound waves. Glasses seem to have a sizable excess of vibrations at this boson peak. At present, there is no agreement which forces drive these extra vibrations into the low-frequency region, though several possible explanations have been proposed [1][2][3][4][5][6][7][8]. Another controversial question [9,10] is whether the interaction with these vibrations is the physical reason for the crossing of the Ioffe-Regel limit, the reduction of the mean free path down to below the wavelength of sound waves in the THz range.

The present paper proposes a detailed atomic model for the boson peak modes in glasses consisting of close-packed atoms or spherical molecules: The soft modes are ascribed to small regions with a pronounced shear misfit with respect to the surrounding matrix. In the following section II we describe the picture and derive the properties of these modes. In section III we relate our model to experiment and in section IV discuss and summarize the results.

II. THE GLIDING-TRIANGLE MECHANISM A. Shear strain defects

The central concept of this model is a structural rearrangement of a limited region in the sample which changes its shape to a sheared one. In three-dimensional close packing, the lowest-energy structural rearrangement can be viewed as a gliding of a triangle of three close-packed atoms over an underlying close-packed plane (Fig. 1). On the left side of Fig. 1, the center of the upper triangle lies exactly over the center of a lower triangle within the close-packed plane. The central points of the six atoms form an octahedron. Gliding the triangle in vertical direction one reaches again a stable position after a distance of 2r/ √ 3 (r sphere radius), right part of Fig. 1. The six atoms there form a pair of edge-sharing tetrahedra. Since the distance between close-packed planes is 2r 2/3, the shear angle is 1/ √ 2 in radians, i. e. 40.5 degrees.

The physical problem of a small piece of matter able to transform to a sheared shape D. Eshelby [11,[START_REF] Mura | Micromechanics of Defects in Solids[END_REF]. Here, we translate Eshelby's result into the usual convention, in which the shear angle e and the shear stress σ are related by σ = Ge (G infinite frequency shear modulus). Let v i be the volume of the spherical inclusion and e i the shear angle difference between its two stable configurations (in the example of Fig. 1, e i = 1/ √ 2 and v i = 6v, where v is the atomic volume). For strictly harmonic potentials, the symmetric case of two equally strained configurations has the elastic deformation energy

E a = γ 8 Gv i e 2 i . (1)
The coefficient γ is given by

γ = 7 -5σ P 15(1 -σ P ) , (2)
where σ P is Poisson's ratio. Since Poisson's ratio lies between 0.1 and 0.44 for the known glasses, γ lies between 0.48 and 0.57, close to 1/2.

Eshelby's solution divides the energy into two almost equal parts, one located in the inclusion and one outside. Their ratio is γ/(1 -γ). In the symmetric case, the inclusion would have to distort by e i /2 to fit exactly into the unstrained hole in the surrounding matrix. From a structural point of view, this is the saddle point between the two stable configurations.

In strictly harmonic approximation, Eq. ( 1) holds and the structural saddle point indeed is a saddle point of the energy landscape. However, the harmonic approximation does not account for the flattening at the saddle point, at least as regards the distortion energy E c of the inclusion. For the inclusion, a better description is by a cosine law

E c = 3Gv 4π 2 (1 + cos 2 √ 2πe c ) (3)
shown in Fig. 2, where e c is the shear angle of the central inclusion, counted from the saddle point of the energy of the inclusion. This cosine law has the correct second derivative of 6Gv at e c = ±1/2 √ 2, the two structural energy minima of Fig. 1.

The embedding energy outside the inclusion depends on the difference between e c and the zero point e m of the shear of the matrix. In terms of e c , the shear state of the hole after removing the inclusion and relaxing the matrix is e c = e m . Let us assume γ = 1/2. Then the second derivative of this outside embedding term must also be 6Gv.

E el = 3Gv 4π 2 (1 + cos 2 √ 2πe c ) + 3Gv(e m -e c ) 2 . ( 4
)
Note that the symmetric case at e m = 0 (the red lines in Fig. 2) has its minimum at e c = 0, with all the energy inside the inclusion and a vanishing restoring force constant for a displacement of e c . Thus one does not get a saddle point of the energy landscape, but one gets a shear mode with zero restoring force constant at the structural saddle point. In this case, there is no elastic energy outside the inclusion and the creation energy E s for the soft mode configuration is

E s = 3Gv 2π 2 . ( 5
)
In order to know the probability of finding such a soft gliding-triangle mode in a closed packed glass, one needs to know the ratio Gv/k B T g , where T g is the glass transition temperature, at which the glass falls out of the thermal equilibrium. For metallic glasses, this information can be extracted from a recent data collection [START_REF] Johnson | [END_REF]. Fig. 3 plots the product Gv for 30 metallic glasses at the glass transition as a function of T g . There is a marked scatter, but on average Gv/k B T g = 17.6 (the line in Fig. 3). With eq. ( 5), this implies a ratio

E s /k B T g of 2.
67, i.e. a Boltzmann factor of 0.069. Thus one has to reckon with a sizeable number of nearly unstable resonant shear modes in a metallic glass.

B. Properties of soft shear modes

The equilibrium value e c0 of the inner coordinate e c for a fixed e m0 is given by the zero of the first derivative of E el , Eq. ( 4), with respect to e c :

e m0 = e c0 - sin 2 √ 2πe c0 2 √ 2π . (6)
The second derivative of Eq. ( 4) at e c0 gives the restoring force constant

D e = ∂ 2 E s ∂e 2 c = 6Gv(1 -cos 2 √ 2πe c0 ), (7)
which vanishes for e m0 = e c0 = 0.

To calculate the kinetic energy of the mode, consider a sphere with volume of 6v, corresponding to the six atoms in Fig. 1. Its radius r 0 is given by

4πr 3 0 3 = 6v. (8) 4
where ρ is the density. The kinetic energy outside the sphere, in which the r 4 of Eq. ( 9) is replaced by r 6 0 /r 2 and r is taken from r 0 to infinity, is a factor of five higher than the one inside. Replacing the density ρ = M/v, where M is the atomic mass, one finally arrives at

E kin = 1 2 2 1/3 3 10/3 5π 2/3 M v 2/3 ėc 2 . ( 10
)
For many purposes, it is useful to define a normal coordinate, A, of the soft mode with

E kin = Ȧ2 /2.
For the soft shear modes

A = 2 1/6 3 5/3 5 1/2 π 1/3 M 1/2 v 1/3 e c . (11)
Consider the shear mode at the structural energy minima, say the octahedron of Fig. 1 at e m = -1/2 √ 2. Then, D e = 12Gv, and the mode frequency attains its maximum ω max with 2 1/3 3 10/3

5π 2/3 M v 2/3 ω 2 max = 12Gv. ( 12
)
It is interesting to compare this frequency with the Debye frequency ω D , determined essentially by the transverse sound modes. The transverse sound velocity v t is given by

v 2 t = G ρ = Gv M . (13)
The Debye sound velocity v D is given by the average

3 v 3 D = 1 v 3 l + 2 v 3 t , (14)
where v l is the longitudinal sound velocity, on the average for the metallic glasses [START_REF] Johnson | [END_REF] a factor of 2.26 higher than the transverse one. For this case, v D = 1.13v t . The Debye wave vector is k D = (6π 2 /v) 1/3 , so on average

ω D = 4.4 v t v 1/3 = 2.72 ω max . (15)
This shows that the gliding-triangle mode considered here has a maximum frequency of less than half the Debye frequency in the absence of any structural strain. In order to estimate the density of soft shear modes, one first needs the number of octahedra per atom in the disordered structure. In the crystalline close-packed structures, both fcc and hcp, one has one octahedron per atom. Though this number might easily be smaller in the disordered case, we will take it there also as one.

For each octahedron, there are twelve possibilities to distort to a pair of edge-sharing tetrahedra. There are four choices of the two triangles which are to glide against each other and for each pair there are three glide directions, 120 • rotated to each other. Consequently, there are two stable minima per atom, the octahedron and the edge-sharing pair of tetrahedra, and twelve saddle points between them. As shown in the Appendix, the twelve saddle points lie pairwise rather close to each other, with an angle of only 38.9 degrees between the two members of each pair.

At the glass temperature T g , the saddle points are already sufficiently high in energy to be, in first order, negligible for the partition function, which can be approximated by a harmonic potential around the octahedron energy minimum at e oct = -1/2 √ 2. Close to the minimum, Eq. ( 6) gives e c0 -e oct = (e m0 -e oct )/2. At the energy minimum, Eq. ( 4) takes the simple quadratic form

E el = 3Gv 2 (e m0 -e oct ) 2 , (16)
with the partition function

Z = 2πk B T g 3GV (17)
at the glass temperature T g .

The probability p s = p(e m = 0) to find e m close to the structural saddle point is

p s ≈ 3 exp(-E s /k B T g ) Z = 54Gv πk B T g exp -3Gv 2π 2 k B T . ( 18
)
The factor 3 stems from the fact that there are three saddle point pairs per minimum, as explained in detail in the Appendix.

It is useful to introduce the coordinate

x x = √ 2πe c0 (19)
which is zero at the structural saddle point and ±π/2 at the two energy minima. From equs.

As one moves away in e m from the structural saddle point, the probability density changes only slowly, but x and the restoring force constant of the mode change rapidly (see Fig. 2).

The probability density of x is related to the one of e m by the derivative ∂e m /∂x calculated from eq. ( 6)

∂e m ∂x = √ 2 π sin 2 x. (21)
The vibrational density of states of the soft shear modes is given by

g • s (ω) = p s 3 ∂e m0 ∂x ∂x ∂ω = p s ω 2 3 √ 2πω 3 m 1 1 -ω 2 /ω 2 m ( 22
)
where the prefactor 1/3 results from the usual convention of normalizing the three modes per atom to 1. Since we are only interested in the density of states well below ω m , we neglect the square root term. Then g s (ω)/ω 2 is simply constant

g • s (ω) ω 2 = p s 3 √ 2πω 3 m . ( 23
)
This constant is to be compared with the constant 3/ω 3 D , obtained by dividing the Debye density of states by ω 2 . For the average Gv/k B T g = 17.6 (see Fig. 2), one finds with Eqs. ( 15) and ( 18) a ratio of 0.6 between the density of the soft shear modes and the Debye one. Thus one finds a rather large number of soft gliding triangle modes at low frequency, comparable with the number of sound waves.

D. The shape of the boson peak

The soft shear modes are not exact eigenmodes of the harmonic vibrations of the glass.

They interact with the other vibrations and in particular with the sound waves. The low frequency shear modes can be understood as the cores of quasi-localized or resonant vibrations. They are bilinearly coupled to the other modes. This system of localized soft vibrations coupled bilinearly to extended modes (sound waves) is the basis of the soft potential model [14][15][16].

Through their interaction with the sound waves, there is also an interaction between the soft shear modes. The bilinear interaction of the soft modes with the much larger number of higher frequency modes causes a downshift of the soft modes, some of which even become harmonically unstable. Stabilization via anharmonicity results in a soft mode spectrum linear in ω below a frequency ω c given by the interaction strength, while at low frequency one has a DOS ∝ ω 4 . The crossover between the ω 4 and ω regions then leads to the boson peak. [17] The total number of soft modes is not changed by these frequency shifts. It has been shown that this interaction mechanism can explain the nearly universal density of observed two level systems in glasses. [18] For the purpose of the present paper, it suffices to know that the exact shape of the boson peak requires a more detailed theoretical treatment. Here, we restrict ourselves to the conclusions which are easily accessible from our postulate for the eigenvector of the excess modes. As will be seen in the next subsection, these conclusions include a calculation of the damping and softening of the sound waves from the measured spectrum of a given metallic glass.

E. Coupling, damping and softening

The coupling of a gliding-triangle soft shear mode to an external shear strain is given by the derivative of the elastic energy of eq. ( 4) with respect to both the external strain e m and the internal coordinate

e c ∂ 2 E el ∂e m ∂e c = -6Gv (24)
if the external shear strain is in the direction of e c . In the general case, the directional average over the five possible shear strain orientations must be taken.

It is usual to express the coupling in terms of a product Λ t Ae m of the shear strain and the soft mode coordinate A, defined via eq. ( 11). From Λ t , one can define a frequency ω t via

ω 2 t = Λ 2 t M v 2 t . (25)
Sorting all the factors out, one finally gets

ω t = 2 5/6 π 1/3 3 2/3 v t v 1/3 , (26)
about a factor of 3.5 lower than the Debye frequency. With ω t , the damping Γ of the transverse sound waves (the full width at half maximum in a Brillouin experiment) is given by

Γ t = 3π 2 ω 2 t g s (ω), (27)
so the treatment predicts a damping of the transverse sound which increases with the frequency or wavevector squared. The Ioffe-Regel condition Γ = ω/π is reached at the frequency

ω JRt ω JRt = 2ω 3 D 9π 2 αω 2 t . ( 28
)
Here α is the ratio between Debye density of states and soft shear wave density of states, which in the average metallic glass should be close to 1. For the average case with Gv/k B T g = 17.6, ω JRt lies at 0.23 ω D .

The Ioffe-Regel limit for the longitudinal sound waves lies markedly higher, because the gliding-triangle modes couple only to the shear (or essentially so; we come back to this point in the discussion). Since the longitudinal elastic constant is given by

C 11 = B + 4G/3,
where B is the bulk modulus, a pure shear coupling implies Λ 2 l = 4Λ 2 t /3 for the longitudinal coupling constant. With again ω 2 l = Λ 2 l /M v 2 l , ω 2 l will be a factor of 4v 2 t /3v 2 l smaller than ω 2 t , with v l /v t ≈ 2.26 a factor of 0.26 smaller [START_REF] Johnson | [END_REF], which makes ω JRl a factor of nearly 4 higher than ω JRt .

The influence of the soft shear modes on the shear modulus is given by

∆G G = 3ω 2 t ωmax 0 g s (ω) ω 2 dω, ( 29
)
where ∆G is the difference between G above and below the boson peak.

III. COMPARISON TO EXPERIMENT

A. Vit-4, a heavily studied metallic glass

The best candidate for a check of these ideas is the metallic glass Vit-4, Zr 46.8 Ti 8.2 Cu 7.5 Ni 10 Be 27.5 , a glass with an exceptionally low critical cooling rate [21], for which every relevant quantity has been measured [19,20,22,23], though the mechanical [23], density [22] and neutron [19,20] data were all taken for slightly different compositions.

At the glass temperature T g = 615 K, the shear modulus is 35 GPa, the bulk modulus 112 GPa [23] and the atomic volume is 0.017 nm 3 , so ω D = 27.5 meV. Fig. 4 shows the boson peak at room temperature determined from inelastic neutron scattering data [19,20] in units of the Debye density of states (at room temperature [23], ω D = 28.6 meV).

If one takes the excess over the Debye density of states as gliding-triangle modes, one can calculate the corresponding damping of the sound waves. One finds that the Ioffe-Regel limit is never reached, even for the transverse sound waves. At the boson peak, the damping is about half the Ioffe-Regel limit. Since the longitudinal waves are even less damped, one expects rather sharp longitudinal Brillouin peaks in an inelastic x-ray experiment. This was in fact measured for another metallic glass [10], Ni 33 Zr 67 . In this case, one finds a Q 2 dependence above ω max , which probably has nothing to do with the boson peak.

Though the boson peak is small, it is crucial for the mechanical properties and their temperature dependence in the undercooled liquid. One calculates ∆G/G = 0.08 from eq. ( 29) with the excess modes from Fig. 4. This fits in with the measured difference [23] of 26 % of the shear moduli of glass and crystal. The shear modulus Grüneisen parameter Γ G = ∂ ln G/∂ ln v is found to be 6.5 from the temperature dependence of G [23] and the density [22] in the glass phase. Consequently, the density ratio 1.027 between crystal and glass explains two thirds of the shear modulus difference, the gliding-triangle modes explain the rest.

In the undercooled liquid, the number of gliding-triangle modes increases exponentially according to eq. ( 18) (replacing T g by T ). They do not only increase their number, but also shift to lower frequency, as seen in Fig. 4. However, if one makes a rough estimate by a linear extrapolation of the liquid data in Fig. 4 to the frequency zero, the increase of the area under the curve follows in fact eq. (18). Taking ∆G/G = 0.08 as the effect of the gliding-triangle modes frozen in at T g , eq. ( 18) predicts an additional reduction of G of 0.023 GPa/K in the melt. This is almost exactly found in experiment [23] (see Fig. 5). The reduction in the liquid is a factor of 2.4 larger than the Grüneisen expectation on the basis of the thermal expansion [22]. The difference can be frozen into the glassy state by quenching. It is obviously due to configurational changes of the melt. The finding is of great importance for the understanding o f liquids. Our model gives a quantitative physical explanation for it: it is due to shear instabilities, stabilized by the surrounding matrix and requiring a formation energy of about three times k B T g , which reduce the shear modulus.

The experiment [23] reveals also an additional decrease of the bulk modulus with increas- ing temperature, though not quite as pronounced as the one of the shear modulus. If this is also attributed to our shear defects, they must couple to the compression. From the data [23], the coupling constant ratio Λ l /Λ t of these defects should be 1.69 rather than the value 4/3 for a pure shear defect, close to the value 1.6 found for the tunneling states in glasses [24]. We will come back to this point in the next section.

B. Low temperature anomalies

The gliding-triangle model makes some rather detailed predictions for the low temperature glass anomalies [START_REF]Amorphous Solids: Low temperature properties[END_REF]. We show this employing the formalism of the soft potential model, the standard model in the field for which an extensive comparison with experimental data already exists [START_REF] Parshin | [END_REF][27][28].

The soft potential model [START_REF] Parshin | [END_REF][27][28] connects the boson peak modes to the tunneling states [START_REF]Amorphous Solids: Low temperature properties[END_REF] and to low-barrier relaxation. It postulates a constant density of first and second order terms D 1 and D 2 around the value zero for the potential

E sof t = W (D 1 x s + D 2 x 2 s + x 4 s ). (30)
The coordinate x s relates again to the coordinate A of the soft mode with E kin = Ȧ2 /2. Let us denote the fourth order term in this coordinate by v 4 A 4 . Then the zero-point energy of the purely quartic potential is by definition W , given by the equality of kinetic and potential

energy at A 0 W ≡ 2 2A 2 0 = v 4 A 4 0 . (31)
This relation defines both W and A 0 . The kinetic energy results from the uncertainty principle at a confinement within ±A 0 . The energy W separates tunneling states below W from vibrational states above. The coordinate x s = A/A 0 .

The gliding-triangle modes of the present paper are also distributed around a purely quartic potential (the continuous red line in Fig. 2). Its fourth order term 2π 2 Gve 4 c , obtained by the fourth derivation of eq. ( 4), leads in combination with the definition of A via eq. ( 11) to W = 5 2/3 π 10/9 2 5/9 3 20/9 (Gv)

1/3 ( 2 /M v 2/3 ) 2/3 . (32)
Eq. ( 32 with the soft potential model for a number of nonmetallic glasses [27]. It combines a very high energy, Gv, coming from the fourth order potential term, with the very low kinetic confinement energy 2 /M v 2/3 . In this context, it is interesting to note that the fourth order term of the quartic potential in Fig. 2 does not result from the short range atomic repulsion, but rather from the decrease of the curvature of the inclusion energy with increasing e c , an anharmonic property of the saddle point of the inclusion energy.

Let us next turn to the tunneling states. Of course, the potential of eq. ( 4) is in no case a double-well potential with a low barrier able to form a tunneling state, because the negative force constant of the inclusion can at most cancel the positive force constant of the surrounding matrix. In terms of the soft-potential model, the treatment in Section II corresponds to the case where the density p(D 2 ) is not constant around zero, but rather a δ-function at zero (this explains why one gets an ω 2 -rise in the excess density of states rather than the ω 4 -rise of the soft potential model). But it was already pointed out in section II that the interaction between the quasilocalized modes is responsible for the shape of the boson peak and the appearance of tunneling states [17,18]. However, such a detailed modelling is beyond the scope of the present paper. Here we make no statement on the number of these tunneling states. But we can make a statement on their coupling to the sound waves.

For W = 4.4 K, a tunneling state with a splitting around 1 K requires D 2 = -6. The two minima of the double-well potential have a distance √ -2D 2 in the coordinate x s of eq. ( 30)

from each other. Using the definitions of equs. ( 31) and (11), one gets the corresponding distance in e c . Multiplying this distance with the coupling factor 6Gv of eq. ( 4) and again averaging over the five possible shears, one gets the coupling constant γ t of the tunneling state to a transverse sound wave

γ t = 2 10/9 3 4/9 5 1/3 π 2/9 2 /M v 2/3 Gv 1/6 -D 2 Gv. ( 33
)
For the data collection on metallic glasses [START_REF] Johnson | [END_REF] this yields an average value of 0.49 eV, not too far away from the average value 0.39 of direct measurements of γ t of the so-called lowtemperature "tunneling plateau" in the sound absorption of 18 different inorganic glasses [24]. In this field, scientists have always been wondering what kind of mode coordinate would be able to give such a strong coupling. Here we have for the first time a detailed answer to this question, at least for an important subgroup of glasses. The coupling to longitudinal sound waves is described by the analogous coupling constant γ l . From the consideration at the end of section II on a pure shear defect like the one proposed here, one would expect a ratio γ l /γ t = 4/3 = 1.155. This does not agree with the experimental finding [24] γ l /γ t = 1.6. As pointed out in the previous subsection, the temperature dependence of the elastic moduli in the Vit-4 melt [23] suggests a factor of 1.69, so metallic glasses in this respect do not differ from the general case. At present, we see no obvious reason for this additional coupling to the compression.

IV. DISCUSSION AND SUMMARY

As seen in section III, the Ioffe-Regel limit for longitudinal sound waves for metallic glasses should be by a factor of two to four higher than the one for transverse sound waves.

In fact, an inelastic x-ray Brillouin experiment in a metallic glass reveals a relatively low longitudinal sound wave damping [10].

However, the metallic glass result [10] is an exception rather than the rule [9]. In the nonmetallic glasses, one usually does find the Ioffe-Regel limit close to the boson peak, maybe slightly above the boson peak. How is this possible?

The question is related to three ratios: The ratio of the boson peak height to the Debye level, the ratio of the coupling constants Λ 2 l /Λ 2 t , which for a pure shear defect is 4/3, and the ratio v 2 l /v 2 t . As it turns out [START_REF] Johnson | [END_REF], this last ratio is about 5.1 in the metallic glasses, unusually high (compare silica, where this ratio is 2.3). This and the weakness of the boson peak in many metallic glasses are the two reasons for the weak longitudinal damping. A simulation [29] in a soft sphere glass, where the ratio v 2 l /v 2 t was even higher than for real metallic glasses, showed a very weak damping of the longitudinal waves.

There seems to be some as yet unexplained coupling of the boson peak modes to the compression. As pointed out in section III, the ratio Λ 2 l /Λ 2 t seems to be 1.6 2 = 2.56 rather than 4/3. This feature is missing in our model. If one ascribes it to different packing factors for octahedron and Pair of tetrahedra, Fig. 1, these would have to differ by about 18 %. This explanation seems not to be very convincing. We rather think that the absence of the compression coupling in our model comes from the neglect of the atomic roughness of the shell surrounding the six central atoms, which introduces a random coupling to the compression. Despite this small uncertainty, the evidence for the validity of the gliding triangle model in section III shows that it is able to compete with other explanations of the boson peak [2][3][4][5][6][7][8]. There is no adaptable parameter; one only needs the shear modulus, the atomic mass and the atomic volume. The picture provides a microscopic basis for the empirical soft-potential model [1,[START_REF] Parshin | [END_REF][27][28], not for all glasses, but for the subgroup of metallic glasses.

The gliding triangle motion considered here resembles the "shear transformation zones" invoked for the explanation of the yield stress in metallic glasses [START_REF] Johnson | [END_REF], because both are based on the concept of a Frenkel gliding of adjacent atomic planes for a finite volume. However, the shear transformation zones do not relate to low-frequency vibrations, but to relaxations over a high barrier, the elementary flow process, probably with a much larger core than the six atoms considered here.

In our derivation we have treated the soft shear motion in a mean field approximation for the surrounding material, approximated by an elastic continuum. When the melt is cooled down through the glass transition, correlations between adjacent groups of atoms gain importance and will couple different modes. This can easily be observed in a crystalline model.

Let us take a simple hcp-lattice, with a stacking sequence of ababab... The rearrangement in Fig. 1 then corresponds to a local change of stacking for the three atoms to abacab... A simple simulation shows that this causes two lines of compression along close packed direction emanating from the triangle. The simple triangular shear thus becomes more collective, the motion more correlated. Allowing the symmetry to be broken, as will always be the case in a disordered medium, one of the two compressional chains gains on the expense of the other which will merge to the first one leaving some "free" excess volume to its side. This picture then resembles the chain-like structures of the soft modes found in simulations of soft sphere glasses [30,31] and other metallic systems. These chains resemble extended interstialcies with an atom squeezed into the densely packed direction and excess volume smeared out on its side. This picture is the basis of the intersticialcy model of Granato and coworkers [32].

But even with such a modification, the crucial point of our paper remains unchanged: One needs local instabilities to create a low frequency or tunneling mode in a stable solid.

These instabilities require a strong local shear misfit, which costs energy, an energy of several times the thermal energy at the glass transition. In model descriptions by random matrices or elastic constant fluctuations, e.g. [2,3] the condition of near instability of some vibrational eigenmodes has to be imposed at the same time putting some ad hoc limits to The validity of these two assumptions would explain the universality of the low temperature anomalies in glasses, including those at slightly higher temperatures where the boson peak modes dominate.

V. APPENDIX: THE 12 SADDLE POINTS IN CUBIC NOTATION

In an fcc close packed crystal, the six atoms of the octahedron have the coordinates (±a/2, 0, 0), (0, ±a/2, 0) and (0, 0, ±a/2) in units of the lattice constant a, related to the atomic volume by a 3 = 4v.

Let us first consider the gliding-triangle motion of the two triangles in the (111)-plane along < 211 > in Fig. 6. The upper triangle ABC moves opposite to the lower triangle DEF. This gliding motion involves not only a shear, but also a rotation. To get the rotation-free shear, let us define the two unit vectors a 1 =< 211 > / √ 6 and a 2 =< 111 > / √ 3. With the origin of the position vector r at the center of the octahedron, the displacement vector u( r)

of the corresponding pure shear is given by

2 e c u = ( r a 2 ) a 1 + ( r a 1 ) a 2 , (34)
where e c is the shear angle of the gliding triangle mode of Fig.

E el = (C 11 -C 12 + C 44 )ve 2 c = (2C ′ + C 44 )ve 2 c , (35)
where C ′ = (C 11 -C 12 )/2 is one of the two shear moduli of the cubic crystal and C 44 is the other. Equating both of them to the isotropic shear modulus G of the glass, we are back to the curvature of eq. ( 3) at the octahedral site.

Two thirds of the elastic energy of eq. ( 35) come from the diagonal elements, a dilatation in the y and in the z-direction accompanied by a factor of two larger contraction in the In the five-dimensional shear space, the lines connecting those two saddle points with the octahedron have an angle of only 38.9 degrees (Fig. 7).

The closeness has consequences for the probability density p s of section II. This was calculated integrating the partition function Z only over the connection line between octahedron and saddle point, and then multiplying with the number of saddle points. In principle, one needs the full integration over the five-dimensional shear space. As long as the saddle points are well separated in phase space, the integration over the irrelevant degrees of freedom can be omitted [START_REF] Vineyard | [END_REF]. But if two saddle points are as close together as in Fig. 7, the integration over their neighborhood (along the dashed lines in Fig. 7) leads to double-counting. In a crude approximation, we counted the four saddle points with the common diagonal shear elements as only two in eq. ( 18) (remember that the equation is only needed to demonstrate 16 determined from inelastic neutron scattering [19,20] in units of the glass Debye-level 3/E 3 D . The arrow marks the maximum frequency ω max up to which one expects gliding triangle modes. Fig. 5: Shear modulus measurement [23] in glassy and liquid Vit-4. The dotted line corresponds to a Grüneisen parameter Γ G = 6.5 in the glass.

  matrix has been treated fifty years ago by several authors, notably by J.

sin 3

 3 Inside the sphere, the motion is a pure shear with amplitude e = e c -e c0 , say u x = ey/2, u y = ex/2 and u z = 0. The amplitude of the outside motion decays with the squared distance from the center. It joins the inner shear amplitude continuously at the sphere boundary. The kinetic energy inside the sphere is E kin,i = ρ ėc θ dθdφr 4 dr,

  of gliding-triangle modes

( 7 )

 7 and[START_REF] Mura | Micromechanics of Defects in Solids[END_REF], one obtains for the frequency of the shear mode ω = ω max sin x.

  ) supplies an average value W/k B of 4.4 K for the thirty metallic glasses of the data collection of Johnson and Samwer[START_REF] Johnson | [END_REF], a value well within the range of those fitted 11

12

 12

13

 13

14

 14

  actual instability. The physical origin of both the disorder and the afore mentioned conditions imposed on the randomness remain open in these models. Our model, on the contrary, deals with the underlying physics -how an inherent structural softness leads to a boson peak.The question is: To which extent is this answer universal? How should one generalize it to other glasses, glasses with covalent bonds forming a random network, molecular glasses with the additional degree of freedom of the molecular rotation, to the random chain structure of polymers? The answer to the question requires further studies. It might be an intelligent guess, however, to assume the following two general features:(i) The boson peak and the tunneling modes are always of the same nature, as postulated by the soft potential model (ii) Boson peak modes and tunneling modes are due to small regions with a strong elastic shear misfit, lying close to a saddle point of the inclusion energy between two stable structural minima of the inclusion.

  x-direction. As it turns out, this special shear is a common feature of four of the twelve saddle point distortions, each coming from a different triangle plane. They only differ in the nondiagonal elements, the one from the (111)-plane having e xy = 1, e xz = -1 and e yz = -2, the one from the (11 1)-plane having e xy = -1, e xz = 1 and e yz = -2 and finally the one from the ( 111)-plane having e xy = 1, e xz = 1 and e yz = 2. This last one is closest to the first one, because the shear occurs in the same shear plane, with only a slight rotation of the shear directions.

Fig. 1 :Fig. 2 :

 12 Fig. 1 :The gliding-triangle rearrangement of six closely packed spherical atoms. Fig. 2: Potential energies for inclusion and matrix as a function of the inclusion shear e c around the structural saddle point, where the matrix equilibrium shear e m lies in the middle between the octahedron and the two tetrahedra of Fig. 1. The continuous black line is the inclusion energy, the dashed red line and the dot-dot-dashed blue line show the matrix energy for e m = 0 and 0.05, respectively. The red dot-dashed and blue dotted lines show the corresponding total energies.

Fig. 3 :

 3 Fig.3: Energy Gv as a function of the glass transition temperature T g in metallic glasses[START_REF] Johnson | [END_REF]. The line depicts the average Gv/k B T g = 17.6.

Fig. 4 :

 4 Fig. 4: The boson peak of the model glass Vit-4 in comparison to crystal and liquid,

Fig. 6 :

 6 Fig.6: The octahedron and one of its twelve gliding triangle modes in cubic notation.Triangle ABC glides in the < 211 > direction, triangle DEF in the opposite direction.

Fig. 7 : 19 Page

 719 Fig. 7: The two neighboring saddle points with a strong common shear e common and a weak opposite shear e xy + e xz with respect to the undistorted octahedron.

  Fig. 1 86x57mm (600 x 600 DPI)

  Fig. 2 101x66mm (600 x 600 DPI)

  Fig. 3 97x66mm (600 x 600 DPI)

  Fig. 4 88x70mm (600 x 600 DPI)

  Fig. 5 93x68mm (600 x 600 DPI)

  Fig. 6 79x57mm (600 x 600 DPI)

  Fig. 7 94x68mm (600 x 600 DPI)

  From eq.[START_REF] Vineyard | [END_REF], one can calculate the shear. In Kittel's notation[START_REF] Kittel | Introduction to Solid State Physics[END_REF] (diagonal elements e xx = ∂u x /∂x, nondiagonal elements e xy = ∂u x /∂y + ∂u y /∂x), one finds for the gliding triangle mode of Fig.6in units of e c √ 2/6 the diagonal elements e xx = -4, e yy = 2 and e zz = 2, the nondiagonal elements e xy = -1, e xz = -1 and e yz = 2. The elastic energy of the octahedron is

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y
	6.
	15

Page 16 of 26 http

  that the density of soft gliding triangle modes is comparable to the one of the sound waves).

	Philosophical Magazine & Philosophical Magazine Letters
	F o r F o r F o r
	P P P
	e e e
	e r e r e r
	R R R
	e e e
	v i e v i e v i e
	w w w
	O n l O n l O n l
	y y y

://mc.manuscriptcentral.com/pm-pml

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters


Maryse Lancin
email: marsye.lancin@univ-cezanne.fr

Michael Texier
email: michael.texier@univ-cezanne.fr

Gabrielle Regula
email: gabrielle.regula@univ-cezanne.fr

Bernard Pichaud
email: bernard.pichaud@univ-cezanne.fr

M Texier

Michaël Texier

Defects created in N-doped 4H-SiC in the brittle regime: Stacking fault multiplicity and dislocation cores

Keywords: CTEM, defect structures, dislocations, FIB, HRTEM, microstructure, quantum wells, SiC stacking fault

    Defects created in N-doped 4H-SiC by flexion in the brittle regime: Stacking fault multiplicity and dislocation cores.

Introduction

SiC was first developed for thermo-mechanical applications because of its high strength, its great wear resistance and its low thermal expansion coefficient. The potentialities of this ceramic semiconductor which possesses a large band gap were not ignored but its applications require a drastic control of the polytype homogeneity, doping level and crystalline perfection of materials. The challenge, now overcome, consisted in growing ingots made up of a sole kind of polytype and containing a few micro-pipes, a dislocation density lower than 10 4 cm×cm -3 and controlled impurity content. With such single crystals, it is possible to carry out fundamental research in order to correlate the dynamics and structural characteristics of the dislocations. Indeed, the dislocation dynamics is only known thanks to controversial theoretical and experimental results. Many experiments on the more common polytypes -3C, 4H, 6H-give the lowest activation energy for silicon core partial dislocations (PD) [1][2][3][4][5][6] whereas other studies show similar mobility for both carbon and silicon core PDs [7]. Note that the leading PDs are identified using the same technique (weak beam-dark field (WB-DF) imaging). In addition, first principle calculations give lowest activation energies of migration for the carbon core PDs [8,9]. Another controversial issue has recently focused interest. It concerns the multiplicity of the stacking faults (SFs) created in the brittle regime in highly N-doped 4H-SiC. Indeed, the identified defects consist of single stacking faults (SSFs) [2,4,5], double stacking faults (DSFs) [7] or multiple stacking faults (MSFs) [7,10]. The observed differences are particularly a matter of concern when the stacking is determined by HRTEM. To account for these discrepancies, different influences of the possible driving forces may be set forth. Three driving forcesthermodynamic, electronic and mechanical -may indeed play a role in these ranges of temperature and doping level (N= 5×10 18 cm -3 -9×10 19 cm -3 ) but their respective weight may vary depending on the deformation process or the N content. It thus appeared interesting to check the role of these driving forces while using a method which allows characterising unambiguously individual stacking faults and their bounding partial dislocations (PDs). To begin with, this work is carried out in N-doped 4H-SiC (2 ×10 18 cm -3 ), which is the most promising material for power electronics and hence was commercially available in wafers of high crystalline perfection. The temperature range of the deformation was chosen according to the literature to develop preferentially Shockley PDs leaving SFs in their wakes. Our knowhow regarding one deformation procedure suitable for semiconductors such as silicon or gallium arsenide was adapted to silicon carbide. Using this procedure, the nucleation occurs at or next to the surface and the defects expand in a well controlled stress field [11]. The density of created defects being low enough, it was possible to characterize individual stacking faults and their bounding PDs by coupling chemical etching, TEM techniques -WB-DF and HRTEM imaging, -with image contrast analysis and dislocation core reconstructions. In a first study performed at 823 K, we have characterised DSFs dragged by Si(g) 3 PDs [12,13]. However, because of the above mentioned discrepancies, we found interesting to perform further investigations paying attention to obtain good statistics on the fault stacking and core composition of the PDs. Moreover, we have included large angle convergent beam electron diffraction (LACBED) to our TEM studies [14] which is a powerful tool to characterise the PD core when coupled with core reconstructions. Besides, we have characterised PDs on each side of the SF intersections with the sample surface and investigated the orientation of the leading PDs versus the applied stress. This paper reports the analysis of all results which allows to discuss the PD mobility versus their core and to understand the features of the DSF expansions.

Materials and techniques

N-doped 4H-SiC was purchased from Cree Research™. The N content was controlled by photoluminescence measurements and found to be 2×10 18 cm -3 [15]. The density of dislocations in the basal plane measured by X-Ray-Topography was lower than 10 4 cm.cm -3 [11].

The 20×5×0.2 mm 3 samples were cut by Cree Research so that their length (X axis) was at 45° from the glide planes. The dislocation sources were obtained by a scratch on the ( )

0 2 11
surface and the defect expansion by annealing under cantilever bending in the brittle regime at 823 K or 973 K [11]. After annealing, the samples were chemically etched and imaged with an optical microscope.

(0001) thin foils for WB-DF TEM and LACBED investigations were prepared by the Focused Ion Beam technique (FIB) at the fault tips. Before ion thinning, platinum was deposited on the sample surface along the fault emergence to avoid any ion erosion. The dislocation lines were imaged by WB-DF TEM with 200 keV microscopes (FEI™ and 200CX Jeol™). The line, direction and modulus of the Burgers vectors were determined by LACBED using a 2010-FEG Jeol™ microscope operating at 197 keV [14].

( )

0 2 11
thin foils for HRTEM imaging were prepared by back-side mechanical polishing and ion-beam milling to electron transparency. They were cleaned with Ar plasma to remove the contamination introduced by the sample preparation. HRTEM studies were carried out with a 2010-FEG Jeol™ microscope equipped with a GIF and a CCD camera (Gatan™) or a 4000EX Jeol™ microscope equipped with a CCD camera (Gatan™). A few focal series were performed [12] but most of the defects were imaged when the minimum of contrast was obtained on the amorphous layer created by the thinning process. Indeed, such a defocus could unambiguously be selected during the experiments.

The simulation of the image contrast was performed using the multi-slice method implemented with the EMS software [16]. The absorption and Debye-Waller coefficients were equal to 0.04 and 0.004 nm 2 for Si and to 0.03 and 0.004 nm 2 for C, respectively. The selected sample thickness (3nm) provided the best fit between experimental and simulated contrasts but it is known that this value is much lower than the true one [17].

The respective positions of Si and C atomic columns in the sample were verified by analysis of the { 0 2 11 } image contrast [18] (figure 1). At the minimum of contrast (resp. Scherzer), the projections in a { 0 2 11 } plane of the Si and C atomic columns which are 0.109 nm apart and form the dumbbells could not be resolved on the images and appeared as white (resp. black) dots. However, they could be localized using the respective positions of the dots in the two sets of glide planes labelled G1G2 and G3G4 (quoted ++ and --respectively in the Hägg notation, with which we think that our results would be less clearly stated). Inside the thick parts of experimental images, G1G2 could be easily distinguished from G3G4 using their "black and white contrast" [19].

Results

Identification of the faults expanding during annealing

The etch patterns obtained on the ( )

0 2 11
face in all the samples annealed at 823 K [12] or at 973 K (figure 2) consist of lines that stretch at 45° from the scratch and thus are parallel to the basal planes. They correspond to the intersections of planar defects with the surface. As previously reported for 823 K [12], three populations of defects labelled A, B and C can be categorised thanks to their propagation direction and their length. The A defects expand as a function of the bending stress in the P1= ] 100 1 [ direction.

They are nucleated next to the scratch and also next to the sample border labelled "edge1" in figure 2.

Notice that the sample edges were not smoothed after the sample cutting and thus could act as defect sources. The B defects also glide in the P1 direction, from the scratch or from the sample edge 1, but their lengths are almost independent of the bending stress and remain smaller than 200 microns. The C defects develop from the scratch in the P2 direction. Their lengths reach at the most 400 microns and decrease when the bending stress increases. Notice that no defect develops on the ( ) (34) have the b-3C stacking, whatever the temperature and the dislocation source location (scratch or sample edge). In the literature, such defects consisting of six cubic Si-C bi-layers being called DSFs, A, B and C defects will hereafter be named DSFA, DSFB and DSFC. Moreover, we always observe DSFs but never SSF and only one multiple fault consisting of 10 cubic Si-C bi-layers. Coupling optical micrographs and HRTEM images, we find that the defects typically develop by batches which are distant of a few tens of microns and which consist of two or three DSFs a few tens of nanometres apart.

To fully characterise the atomic structure of the DSFs, we have localised the projections of the Si and C atomic columns in the ( ) planes. On the contrary, it would look like the b-3C stacking if it were observed in ( )

20 1 1 , (
)

0 1 1 2 or ( ) 0 1 2 1
planes, due to the P6 3 mc structure of 4H-SiC.

Partial dislocations dragging the faults

The PDs bounding the DSFs were characterised next to the ( )

0 2 11
surface thanks to WB-DF imaging and LACBED analysis. The characterisations were generally performed on the mobile segments which border the defects in their expansion direction, but some were also performed on the immobile segments located near the scratch. Two significant features are shared by the three populations of DSFs. Indeed, the WB-DF images always reveal two SF S dragged by two PDs called D1 and D2 at an equilibrium distance of about 100 nm (figure 4). The extinction criterion demonstrates that D1 and D2 bounding each DSF always have the same Burgers vector and thereby are repulsive. Beyond these two similitudes, each population of DSFs exhibits its own characteristics, whatever the temperature or the dislocation source location (scratch or sample edge). In the following, the direction, L, of the mobile PD is oriented from the ( ) Table 1 summarises the number and characteristics of the PDs identified. When looking back to figure 3, it is clear that the PD pairs driving DSFAs glide on successive G1G2 planes resulting in the a-3C stacking whereas those dragging DSFBs and DSFCs share two successive G3G4 planes giving the b-3C stacking. Note that the two other associations of two adjacent glide planes (G2G3 or G1G4) cannot create a DSF but two layers of 2H-SiC.

Core composition

Knowing a PD line and character as well as its type of glide plane -G1 (G2) or G3 (G4) -it is then possible to obtain its core by reconstruction. To better discuss the dislocation formation, we find it more and which glides in the P1 direction either on G1(G2) or G3(G4) planes (figure 6). The boundaries are rebuilt in the ( )

0 2 11
plane perpendicular to the line in order to obtain the structural units of 90° or 30° PDs, characteristic of elementary semiconductors [20,21] and observed in 3C-SiC [22]. Obviously, the reconstructions show that the PD core depends on its character (30° or 90°) but also on its glide plane and thus is related to the a-3C or b-3C stacking. As an example, 90° PD pairs have a silicon core and result in the a-3C stacking when they glide on G1G2 whereas they have a carbon core and induce the b-3C stacking when they drive the fault on G3G4. Note that the upper reconstruction in figure 6a corresponds to the 90° PD pairs dragging some DSFAs and the lower one in figure 6b, to the 30° PD pairs leading the DSFBs.

The possible Burgers vectors of the different PDs may be determined by drawing the Burgers circuit and using the starting finish/right hand convention. For the 90° PD, we found the b modulus and sign since the vector is in the ( ) Based on the reconstructions, all the possible loops which result in the a-3C or b-3C stacking can be sketched (figure 7). As it is well known, in such theoretical loops, the partial segments follow the Peierls valleys. In a compound semiconductor such as SiC, only six loops may expand in the material. Moreover, because the PD sources are located on or very near the ( )

0 2 11
surface, only half loops could develop during the deformation process. One DSF consists of two half loops expanding in successive G1G2 or G3G4 glide planes. In figure 7, the three observed populations of DSFs are indicated.

We think important to underline here the potentialities of the method used to characterise the defects.

The reconstructions allow the PD core determination, knowing the stacking and the propagation direction of the fault as well as the PD characteristics obtained by WB-DF TEM (line and Burgers vector direction).

In addition, and very interestingly, these reconstructions give the PD core and the stacking of the fault when the PD Burgers vector modulus and sign are determined by LACBED. Figure 7 shows that a PD pair

Discussion

The deformation procedure selected induces the nucleation of PD pairs, gliding on adjacent glide planes in which the Si-C dumbbells have the same orientations, resulting in the formation of DSFs, whatever the temperature, the source location or the applied stress. Surprisingly, it promotes half loop glide on only one side of the scratch, either in the P1 or P2 direction. Moreover, it creates only three of the six populations expected according to the reconstructions (figure 7). It is noteworthy that all the identified PDs dragging the DSFs have a silicon core. Most of the mobile PDs are of the 30° type but, interestingly, we found some DSFAs driven by Si(g) segments with a 90° character.

To understand such a deformation process, let us first consider the Peach-Koehler force F applied on each half-loop by cantilever bending during our experiments [11]. We remind that

F = [σ r ] b ^ L (1)
and that the resolved shear stress is given by:

         - = 0 0 0 0 0 0 0 0 k r σ (2)
with k(X,Z) being positive above the neutral plane (zone under compression); k depends on the SiC Young modulus (E), the sample thickness (e), the Schmid factor (s) and the curvature radius R(X). Therefore, three of the possible half-loops (HL1, HL5 and HL6) expand under a compressive stress and the other three (HL2, HL3 and HL4) under a tensile one. With the ( )

0 2 11
face submitted to a compressive stress, we would expect symmetric expansions of HL1, HL5 and HL6 but we only observe an asymmetric glide of the HL1 pairs (DSFAs), HL6 pairs (DSFBs) and no HL5. This implies that the Si(g) segments are much more mobile than the C(g) ones. Under such a condition, a loop containing only one Si(g) segment such as HL5 would hardly grow and cannot be detected. The existence of HL4 pairs (DSFCs), which are not expected on this face and which never extend from the sample edge 2, proves that the deformation induced by the scratch provides a glide force expanding them in the P2 direction; this is similar to what a tensile stress would do. This force however, has limited impact as compared to a compressive stress.

Actually, the expansion of DSFCs is very little and decreases as this compressive stress increases. The At this point, a question arises: when a DSF is created, why it is it always dragged by a pair of identical partials? Despite the complexity of the stress field next to the scratch, it is clear that two PDs nucleated on two adjacent glide planes must relieve either compressive or tensile stress. Let us consider the PDs which glide in the P1 direction to relieve a compressive stress. We have seen that PD pairs must glide either in G1G2 or in G3G4 to form a DSF. In G1(G2), only HL1 expands under compressive stress and thus two identical loops must be created. In G3(G4), if HL5 is created, its expansion is blocked by its immobile C(g) segments, as already mentioned, and, hence, the only possibility to obtain a DSF is to nucleate two HL6. The existence of only identical PDs in a pair is perfectly consistent with the greater mobility of Si(g) as compared to C(g). The lack of mobility of C(g) also prevents the expansion of fronts of two PDs -as observed in 3C-SiC at higher temperatures [22] -which would also result in DSFs without the contradictory effect of the repulsion of dislocations with identical Burgers vectors. Indeed, such fronts would consist of a combination of 90°Si(g) and 30°C(g) or vice versa.

The mobility of Si(g) and C(g) were reported in heavily N-doped (~ 10 18 -. 10 19 cm -3 ) 4H-SiC. They were obtained for similar temperatures and applied stresses than those we used but thanks to different deformation procedures. During compression and indentation tests [5,2], the faulted loops which might expand to relieve the stresses were always dragged by Si(g). However, both Si(g) or C(g) were identified during compression under confining pressure [7]. Moreover, at a higher temperature (1423 K), near the brittle-ductile transition, Chung et al. [23] reported the expansion of DSFs dragged by Shockley PDs in Ndoped (9×10 19 cm -3 ) 4H-SiC after surface damaging and annealing. By reconstruction of the cores of the PDs displayed in figures 2 and 3 of ref [23], we found that the DSFs are necessarily driven by 30°PDs with different core compositions, namely 30°Si(g) and 30°C(g). In view of these discrepancies, the role of the dislocation nucleation appears to be a relevant parameter.

In our samples, almost all the Si(g) identified near the sample surface are of the 30° type. However, results obtained on the DSFAs prove that the deformation process provides enough energy to create 90°Si(g). The mobility of the Si(g) as a function of their character and of the bending stress will be fully discussed elsewhere [24].

During the deformation, C(g) segments can be created, even if they hardly move. Indeed, as an example, we have shown in figure 4f 4f). Though they look more rugged than the Si ones -since they do have a few pinning points -they do not exhibit the marked ziz-zag morphology very characteristic of the C-core PDs created by indentation in 4H-SiC [2] and in 3C-SiC [1] by Pirouz et al. or observed in PiN diode by Twigg et al. [25]. In their deformed samples, Mussi et al [7] found C(g) core PDs with either a straight or a zigzagged-line. Mussi's and our results show that the line morphology of a PD is not a sufficient criterion to deduce its core composition. Now let us compare our results with the predictions of the theoretical computations performed in SiC.

Blumeneau et al [8] showed that the core energy of the PDs is almost independent of the polytype and that their calculations performed on 3C and 2H polytypes can apply to 4H-SiC. They obtained a core energies smaller for the Si(g) than the C(g) and lower for the 30° Si(g) than the 90°Si(g). They also found that the C(g) are more mobile than the Si(g), the 90°C(g) being the quickest. Savini et al. calculated the migration energies of PDs [9]. They used asymmetric core reconstructions, as Blumeneau et al. did, as well as symmetric ones (SR). The latter were the most mobile and were found to be the most stable in heavily Ndoped materials, consistent with our experiments. The activation energy of migration was found to be the smallest for the SR 90° C(g) (0.61 eV) and almost the same for the SR 30°C(g) and the SR 90°Si(g) (1.57 eV and 1.55 eV, respectively). In our experiments, the possible highest formation energy of C(g) did not prevent their formation. However, because the 30°C(g) can be nucleated by our deformation, if they were as mobile as the 90°Si(g) or quicker, we would have observed symmetric expansions of the half loops on each side of the scratch, instead of the asymmetric ones observed. Hence, there are still large discrepancies between experiments and calculations. However, as suggested by Blumenau et al. [8], these discrepancies may originate in the N content of the samples used in the above mentioned experiments. Indeed, the nitrogen atoms which are always located on the carbon sites may create a barrier to kink migration along the C-core partial dislocations and thus decrease their mobility. Velocity measurements in intrinsic 4H-SiC would be useful to verify this assumption.

The variations in the PD mobility as a function of their core compositions explain the asymmetric expansion of the DSFs but not all the deformation features. The different expansions of DSFAs and DSFBs are a consequence of the compressive force due to bending which is, for a given stress, two times greater on DSFA than on DSFB. Though both DSFAs and DSFBs are nucleated by the scratch or the sample-cutting damage, the primary glide system involving DSFAs is the first to be activated by bending and DSFAs are driven far from their origin. Their expansions must be sufficient to relieve the bending stress and to prevent or limit the activation of the secondary glide system which involves the DSFBs. The most intriguing feature of our deformation experiments is the creation of DSFs dragged by two repulsive PDs. The mechanism proposed for multiple SF formation [10] cannot explain the systematic creation of DSFs. If the stress fields were sufficient to promote one dislocation cross slip, there would not be any reason to limit this process to one single event. A second or several others should be promoted leading to MSFs, in particular for higher resolved shear stresses.

The systematic formation of DSFs and the presence of close repulsive PDs imply that extra driving forces, other than mechanical ones, should play a determinant role in our deformation process. These forces must provide the energy necessary to create and expand the faults in the 4H structure and to counterbalance the repulsive forces between the PDs. The energy costs to form stacking faults of various thicknesses have been evaluated by Iwata et al. in pure 4H-SiC [26,27]. They have calculated the stacking fault energy γ sf to create one SF and a second SF next to the first one, in order to form a DSF. They found 17.4mJ/m 2 for the former and 4.7mJ/m 2 for the latter. Moreover, γ sf of the third and fourth neighbouring faults which would increase the cubic layer thickness are equal to 20.3mJ/m 2 and 5.4mJ/m 2 , respectively.

According to these calculations, the energy cost of a DSF is similar to that of a SSF and lower than the one necessary to create two SSFs. Hence, DSFs should be easily created in pure 4H-SiC but multiple faults should not.

In the literature, two driving forces have been proposed for providing the necessary energy to induce DSF formation and expansion. In one model, the 4H 3C phase transformation provides the energy of formation and glide of PDs, hence expanding the corresponding trailed cubic layer. We suggested this thermodynamic force as a possible extra driving force in our experiments [12,28] because it is known that N-doping stabilizes the 3C structure [29][30][31], particularly in the temperature range on consideration.

However, if N-doping stabilizes the 3C structure, one would expect a negative γ sf in N-doped 4H-SiC and thus the formation of large cubic layers of various thicknesses, as long as nucleation sites for PDs are available. Numerous nucleation sites are certainly created by our deformation procedure, but we observe only one multiple fault out of 157 defects identified, which involves the glide of two PD pairs in two neighbouring G3G4 glide plane pairs. Hence, the thermodynamic force cannot be the only cause of DSF formation in our experiments.

The most popular model to account for SSFs and DSFs is the quantum well action (QWA). This is linked to energy levels in the band gap due to the faults, which are deeper for DSFs than SSFs. This model was first proposed by Miao et al. to explain SSF expansion in PiN diodes [32]. Using first principle calculations, they showed that a SSF induces energy levels in the band gap. In N-doped material, for suitable temperatures, the lowest level is below the Fermi level and the energy gain per electron trapped at one SF provides the driving force to extend it. Skowronski and co-workers proposed the QWA to explain DSF expansion in N-doped 4H-SiC annealed at 1423K with or without surface grinding [23,33]. The is a good candidate to explain DSF formation in our samples because in 4H-SiC, for a doping concentration of N= 2×10 18 cm -3 and a temperature of 823 K or 973 K, the Fermi level is above the lowest energy level of the DSF. In a recent paper, Lambrecht and Miao interpreted the gain in energy to expand a faulted unit area as a driving force applied on the PD(s) bounding the fault [35]. They found this force to be about 100 MPa for a SSF and 850 MPa for a DSF. The latter stress is far larger than the resolved shear stress in our experiments. According to these calculations, SSFs as well as DSFs could form and expand in our materials. As stated above, the nucleation sites are numerous enough in our samples to allow the nucleation of two PDs in two successive glide planes of the same type and thus the expansion of DSFs which reduces the crystal energy more efficiently than SSFs.

Conclusion

This paper deals with characterisation of defects which were introduced in N-doped 4H-SiC to study the PD mobility. We develop a powerful approach to individually characterise the stacking faults and the core composition of their bounding PDs. The approach consists in coupling WB-DF imaging and LACBED performed on thin foils milled by FIB at the fault tips, together with HR imaging of the stacking faults, image contrast analysis and core reconstructions.

Core reconstruction of the PDs dragging the six possible faulted-half loops in a compound semiconductor demonstrates that the stacking of the faults and the PD core composition are linked and depend on the dislocation character but also on the glide plane and direction. Only three kinds of faultedhalf loops are identified in the samples after surface scratching and annealing under a bending stress. These half loops only expand in one direction, on one side of the scratch. Moreover, they always glide in pairs on planes of identical crystalline structures (i.e., the same Si=C dumbbell orientation) and form six Si-C cubic bi-layers called DSFs. In the three populations of DSFs identified, each defect is dragged by a Si-core PD pair, the Burgers vectors being the same for a given pair. All these structural characteristics imply that C-core PDs can be created but are immobile at the temperatures and stresses under consideration, (823 K-973 K and 0-100 MPa, respectively). The DSF formation also proves that extra thermodynamic and electronic driving forces play a role in N-doped SiC submitted to our deformation procedure. . The compressive stress increases with X and decreases with Z from σ(X, Z) on the surface to zero on the neutral plane [11]. The white dots correspond to the projections in the image plane of the Si=C dumbbells. The image contrast (left) is simulated using the multislice method (V = 400 kV; Cs = 1 nm; defocus = 8 nm, slice thickness = 3.1 nm, defocus spread = 5 nm, semiconvergence angle = 1 mrad, objective aperture=16 nm -1 ).

The projections of Si (large dots) and C (small dots) atomic columns superimposed to the a-3C stacking are deduced from the image contrast analysis. They are located by the software [16] on the simulated image. For sake of clarity, only the Si stacking sequences are indicated on the drawing. The projections of the glide planes, labelled G1 to G4, are shown.

  border labelled "edge2".The atomic structure of the defects was obtained from high resolution images on [ ] both annealing temperatures, we observe the double faulted layers displayed in figure3. The ABCB stacking typical of 4H-SiC is transformed into either ABCABC or CBACBΑ, which are characteristic of a cubic structure (3C-SiC). We thus label them a-3C and b-3C. It is noteworthy that all the identified A defects (49) exhibit the a-3C stacking whereas the entire B population(18) and the entire

  the image contrast. The simulated contrast of the a-3C structure and the atomic positions are displayed in figure3. Note that a-3C and b-3C are defined in the , the a-3C stacking would be exactly the same, whether obtained in ( )

. 1 [ 1 [..

 11 plane when the PD glides in the P1 direction, and vice versa when it expands in the P2 direction.We selected PDs bounding DSFAs, which expanded during the annealing over 300 microns or more, to avoid any confusion with DSFBs. WB-DF imaging performed on 27 PDs always show that both D1 and D2 exhibit a Burgers vector, b, equal to ± aof some of them gives the sign of b and thus b = a/3 Concerning the mobile PDs, their lines seem to be local stress dependent.For the smallest bending stresses, the PDs are perfectly aligned along two ] are thus true 30° Shockley PDs. For the highest stresses, segments aligned along the ] 20 1Peierls valley, which are 90° segments, can be identified near the ( ) 4b). As the stress decreases (with X or with Z), the PD lines tend to curve and to lie closer and closer to the have a stronger and stronger 30° character (figure4c). Contrary to PDs dragging the DSFAs, the mobile PD pairs bounding DSFBs exhibit similar morphologies and lines. They are rather straight and their lines remain within 10° close to the ] 20 1 Peierls valley (figure. 4d). By means of WB-DF TEM, we obtain b = ± athe PDs identified being thus of the 30° type. The true b value measured by LACBED is a. In a previous work [12], we concluded that b was ± aThe error was due to an incorrect location of the ( ) 0 2 11 surface in the FIB thin foils. The platinum cap always used in this work to protect the surface prevents such an upside down error. Concerning the mobile PDs bordering the DSFCs, they are generally curved near the sample surface (figure 4e), their angle with the ( ) 60° to 40°. As expected, the PDs tend to align along the [ ] they go deeper into the material. WB-DF TEM investigations of the 18 PD pairs analysed give b = ± aBecause of this b orientation and despite the line variations, we admit that the PDs are mainly of the 30° type. The characteristics of the PDs were also investigated near the DSF sources. The morphology of the immobile segments is very similar to that of the mobile ones, whatever their line directions or the DSF types. Figure 4f displays their typical features in the case of PD lines bordering a DSFC. In this case, they are almost stretched along the ]

  consider all possible loops which result in a-3C and b-3C stackings and determine the core composition of the bounding PDs. To draw them, we choose to consider all possible configurations for a PD pair whose line direction is L=

.

  for DSFA is in perfect agreement with the result obtained separately by LACBED. In the case of 30° PDs however, two possible Burgers vectors correspond to the same reconstruction because it is actually the projection of b in the determined. By coupling the results with those obtained by WB-DF TEM, it is possible to deduce the modulus and sign of b. As an example, lets us consider a PD pair leading a DSFB. The reconstruction results in b= a/Of course, we found the b modulus and sign of DSFBs measured by LACBED.

  -core and necessarily glides in G3G4 inducing the b-3C stacking observed separately by HRTEM for DSFBs. Hence, core reconstructions make it possible to match HRTEM and LACBED without any ambiguity.

  propagation also confirms the different mobility of the PDs as a function of their core compositions.

  that the immobile PD segment bounding one DSFC is almost thus consists of 30° C(g) (cf. left part of HL4 in figure7).Moreover, the DSFB expansion only observed in the P1 direction implies that the 30°C(g) segment of HL6 is created and blocks the glide in the P2 direction. Our work is thus in perfect agreement with that of Pirouz et al., as far as the character and core composition of the mobile PDs are concerned.

  (g) morphology however, there is some discrepancy between our work and what is generally admitted. Let us again consider the C(g) segments observed near the sample surface (

Figure 1 : 1 000

 11 Figure captions

Figure 2 :

 2 Figure 2: Planar defects revealed by chemical etching and optical microscopy on the ( ) 0 2 11 surface.

Figure 3

 3 Figure 3: ] 0 2 11 [ experimental HRTEM images of the 4H structure and of the two kinds of faulted stacking created by the deformation procedure at 823 K or 973 K.

Figure 4 :

 4 Figure 4: Lines of the partial dislocations edging DSFs imaged by weak beam dark field near the (0001) axis. The PD lines are imaged in FIB foils prepared at the dislocation intersections with the surface. 4a to 4c display the typical line orientations of mobile PDs bounding DSFAs which glide either under low (~ 10 MPa), high (~ 80 MPa) or intermediate stress values, respectively. 4d shows the mobile PDs of a DSFB and 4e those of a DSFC. 4f) displays the morphology of immobile PDs edging a DSFC near its source.

Figure 5 :

 5 Figure 5: Experimental LACBED pattern of a partial dislocation pair bounding a DSFB (left) and the schematic representation of the Bragg line splitting produced by the dislocations and stacking fault (right). The Burgers vector of the two PDs was determined by indexing the Bragg lines and analysing the splitting using the g.b = m rule [14]. The dashed arrow indicates the position of the PD pair.

Figure 6 :

 6 Figure 6: Core reconstructions in the ( ) 0 2 11 plane perpendicular to the lines of partial dislocation

Figure 7 :

 7 Figure 7: Possible faulted half loops created by the glide of partial dislocations in G1 (G2) or in G3 (G4) planes. Theoretical half loops, whose segments follow the Peierls valleys, which could expand from a scratch performed on the ( ) 0 2 11 surface. Under a compressive stress, half loops HL1, HL5 and HL6 would expand and half loops HL2, HL3 and HL4 would shrink. Inside each half loop, the arrow shows the corresponding b. The half loop sketches corresponding to DSFA, DSFB and DSFC are indicated.

  Figure 1 99x30mm (600 x 600 DPI)

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

http://mc.manuscriptcentral.com/pm-pml

Acknowledgments

Page 12 of 29 http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters They authors cheerfully thank Dr. J Douin who performed the first WB-DF characterisations and H. Idrissi who performed the sample deformation. The authors greatly thank the "Centre Pluridisciplinaire de Microscopie et de Microanalyse de l'Université Paul Cézanne" where the FEI FIB, 200 keV FEI and 200 keV Jeol 2010F microscopes are available. They gratefully thank Pascale Bayle for allowing access to the Jeol 4000 EX microscope as well as the ONERA-Chatillon and the LPM Poitiers where they performed many WB studies on 200CX keV Jeol microscopes.

Defects

Number analyzed

29 52

Table 1: Structural characteristics of the DSFs and of their bounding PDs.

We give the PD line orientation and the PD character that are the most characteristic of each of the three populations of DSFs (the possible line variations are detailed in the text).