Threading Dislocations Piercing the Free Surface of an Anisotropic Hexagonal Crystal : Review of Theoretical Approaches

Inclined threading dislocations (TDs) piercing the oriented free surface of a crystal are currently observed after growth of oriented thin films on substrates. Up to date the unique way to treat their anisotropic elastic properties nearby the free surface region is to use the integral formalism, which assumes no dislocation core size and needs numerical double integrations. In a first stage of the work, a new and alternative approach to the integral formalism is developed using double Fourier series and the concept of a finite core size, which is often observed in high-resolution transmission electron microscopy. In a second stage, the integral formalism and the Fourier series approaches are applied to the important case of a TD piercing the basal free surface of a hexagonal crystal. For this particular geometry, easy-to-use expressions are derived and compared to a third approach previously known for a plate-like crystal. Finally, the numerical interest and the convergence of these approaches are tested using the basal free surface of the GaN compound, in particular for TDs with Burgers vectors c and (a + c).


Introduction
Epitaxial growth of semiconductor materials offers opportunities in combining and modifying structural and electronic parameters [1][2][3][4][5].Unfortunately, the thin crystalline film deposited on the substrate often contains high densities of dislocations piercing the free surface, the so-called threading dislocations (TDs), which affect adversely device properties.At high temperature, a reduction of the density of some TDs can occur via fusion or annihilation reactions.It is therefore of importance to study their intrinsic elastic properties.The aim of the present work is first to investigate the elastic field of an inclined TD piercing the free surface of any anisotropic, semi-infinite crystal.The inclination angle with respect to the normal is denoted .For isotropic crystals, exact solutions are well known for  < /2 [6,7] and  = /2 [8,9].We present a new, alternative approach to the TD elastic field to that developed previously by Lothe et al. [10] who uses the integral formalism and requires in general a set of double integrations.This new approach is developed in terms of double Fourier series and can take account of a finite dislocation core size, which is not possible with the integral formalism.These two approaches are then applied to the important case of a hexagonal crystal (parameters  and ) limited by a basal plane.The particular geometry of a TD parallel to the c axis enables simplified expressions to be derived.These latter are compared to a third approach developed previously for a plate of finite thickness [11] from stress potentials [12].Since wurtzite gallium nitride (GaN) has attractive optoelectronic applications, such as light emitting diode and laser diodes [1,[3][4][5], this compound is used in the present work to perform numerical applications.

The Fourier Calculation Approach.
A description of the elastic field of inclined dislocations lying in an isotropic plate was recently presented from a Fourier calculation method 2 Advances in Condensed Matter Physics Figure 1: Cartesian frames and symbols attached to (a) an isolated inclined threading dislocation (TD) oriented by //   3 in the plane  2  3 and (b), two infinite sets of periodically distributed TDs oriented along .Each of these TDs has a cylindrical core, the section of which is a small square.[20,21].Below, this method is generalized to treat any anisotropic semi-infinite crystal.For the calculation, we use two Cartesian coordinate systems with the same origin , Figure 1(a).The system  1  2  3 is linked to the free surface and oriented with the axis  3 directed outside matter.The direction  1 is any direction taken in the surface.The second system denoted   1   2   3 , is attached to a TD placed and oriented along   3 .This latter system is obtained from  1  2  3 by a rotation  around  1 .The TD terminates in , is oriented by the unit vector //   3 , and has the Burgers vector b.
The first step of the calculation is to consider a biperiodic series of largely spaced and identical emerging TDs.Its stress field can be written as the sum (,  = 1-3) in which   ∞ is the stress field generated in an infinite crystal and   rel is the unknown relaxation stress field that cancels the forces applied along the plane chosen for the surface.To fulfil (1) the idea is to use partial results obtained in [20] for   ∞ and in [22] for   rel .
Let us first consider the calculation of   ∞ .
In the plane   3 = 0, Figure 1(b), the theoretical dislocation lines form a square pattern of black points.These lines are distant by the period  along   1 and   2 .Each dislocation line is oriented parallel to   3 and is associated with a finite core.According to some high-resolution transmission electron microscopy observations and numerical simulations, for example, [23][24][25], a TD core can have a variety of atomic structures.These latter have certainly an influence on the elastic properties very near the theoretical dislocation line [25].However, beyond a few |b| values and thanks to the St. Venant principle, elasticity can be applied without accounting for true limiting boundary conditions around the core.In the present work, for simplicity, the core section is chosen as a square, limited by an infinite prism with a section  2 , limited by portions of planes parallel to planes   1 =   2 = ±/2.Inside the core the dislocation density is taken as b/ 2 , while outside the prism, the dislocation density is zero.Each prism can therefore be described as an extended core dislocation.Around the axis   3 , when  increases and the core size  decreases, the elastic field tends gradually to that of a classical Volterra dislocation.Along the axes  1 and  2 , the periodicity of the defects is  and /cos .
According to [20], the displacement   ∞ and stress fields   ∞ of the infinite dislocation set can be written in the  1  2  3 frame in the form of double Fourier series.The   ∞ stress field can be expressed as (integers ,  ̸ = 0) in which the preexponential terms   (,) can be calculated explicitly with The terms   (,) depend on the elastic constants   taken in  1  2  3 , on the Burgers vector b, on the inclination angle , and on the geometrical parameters  and .However, since their expressions are much too long to be written extensively, they are not given in the present text.
Let us now consider the calculation of   rel .Since this stress field should also be periodic, it can be expressed formally as a double Fourier series.As shown in [22], the formal forms of the associated displacement field can be written as ( and  ̸ = 0) in which the constants   are three unknown complex constants attached to each harmonic term.They depend on the order (, ) of each harmonic term and on the limiting boundary conditions of the problem.The constants   and   are found from the nontrivial solutions of the equations For  = 1-3, the constants   are complex chosen with a negative complex part.The corresponding constants   are normalized solutions of the same equation.In addition For convergence towards  3 < 0, only the negative complex parts of the constants   are retained.The stress field is finally calculated from derivation of ( 4) and application of the Hooke law: The constants   are The unknown constants   are found from the following condition: the addition of the stress fields ( 2) and ( 7) are such that zero stresses act on the plane For each  value, this expression represents in fact a double Fourier series, which must be identical to zero.The only solution is that all its terms (, ) are zero.This leads, for each term (, ), to a system of three linear equations, the solution of which gives the values of the three complex unknowns   .From the knowledge of the   values associated with each term (, ), the total displacement field is obtained by the addition (u ∞ + u rel ), in which u ∞ is calculated from expression (5) in [20] and u rel from expression (4).The total stress field is finally obtained from the addition of expressions ( 2) and (7).
As a consequence, three parameters are required to evaluate the elastic field around an isolated TD: the period  should much larger than |b|, while the size  could be estimated to a few |b|.Of course, a sufficiently large number of terms (, ) must be retained to obtain a sufficiently precise solution.
For the case of a hexagonal crystal,  = 0 and  3 //c, simplified expressions can be obtained, briefly listed below.The axis  1 can be chosen anywhere in the basal plane.
Constants   are the following: Compact expressions are then obtained for the   constants because the terms  1 and  2 are not involved in the calculations.As a result, for brevity, denoting  3 ( = 1-3,  = 1-3) as   and  3 as   , respectively, the   constants become ) [ Similar expressions are found for the particular cases  = 0 or  = 0.The elastic field described by the above Fourier approach can be compared with that calculated with the integral formalism [10].After a brief presentation of this latter approach, we present simplified expressions for a hexagonal crystal (parameters  and ) that has a basal free surface.

The Integral Formalism Approach.
In the integral formalism, the relaxation stress field of a TD dislocation is calculated for convenience in the frame  chosen in [10] Figure 1(a).The components of the Burgers vector are denoted b = (  ,   ,   ).It is generated by a plane fan-shaped distribution of straight infinitesimal dislocations placed in the plane  3 = 0 of the unlimited crystal.This distribution is located in the plane  = 0 and has the Burgers vector density b() = { 1 (),  2 (),  3 ()}.As in [10],  is the oriented polar angle of an infinitesimal dislocation linked to the plane .The field   rel is obtained from double integrations: the first one is needed to compute the energy matrix () related to each dislocation, while the second one is required to sum the contributions of all the planar set of dislocations.
For the particular orientation  = 0, we obtain compact expressions for the image force  that applies on a TD segment  in the plane  = 0 and for the relaxation stress   rel .The integral (24) given in [10] can be greatly simplified.If  is a distance along the TD taken from the emerging point and  a length element, we obtain the simple expression This expression shows clearly that  ̸ = 0 for mixed dislocations, while  = 0 for pure screw (  = 0) or prismatic edge dislocations (  = 0).This explains the trend, for mixed TDs, to have directions systematically misoriented with respect to the c axis, for example, [2].
When  = 0 the density b() becomes simple to express because the calculation shows that the energy coefficient matrix  attached to an infinitesimal dislocation placed in the plane  = 0 is no more depending on its orientation .This constant matrix can be evaluated from expressions (73, 74) in [26] or other similar procedures indicated in [13,27].In passing, some misprints in [26] are indicated in Appendix A.
The Burgers vector of an infinitesimal dislocation has the following simple form: ) in which  11 ,  22 and  33 are the nonzero elements of the energy coefficient matrix .
The relaxed displacement field u rel in a point defined by its cylindrical coordinates ( 1 ,  1 ,  1 ) is then calculated from the addition of all the contributions of the infinitesimal dislocations with Burgers vector b(), in the plane  = 0.
As an example, for the edge case with b(  , 0, 0), the nonzero components of the stress tensor are in which [10] (see Appendix B) For the screw case with b(0, 0,   ), the stress tensor is The integration domain Δ is indicated in [10] with a crossed integration symbol.Therefore, it is not performed on  exactly.For numerical applications, to avoid zero denominators in ( 13)-( 16), an angular domain Δ running from ( 1 + ) to ( 1 +  − ) should be used with  tending to zero.Our experience is that convergence is ensured with  equal to 0.0001, a value sufficiently small to obtain a good precision on the stresses nearby the TD core.

The Stress Potential Approach.
The displacement field regarding a mixed dislocation with //c( = 0) was already expressed for a plate-like hexagonal crystal in [14] but unfortunately this author did not consider a semi-infinite crystal.In the present work, the solution is found in searching for the limit of the displacement field u rel around the upper emerging point when the plate thickness tends to infinity.The calculation is carried out in the  frame, which gives the following result: (i) Edge dislocation b//: the field u rel is derived from two stress potentials  1 and  2 mentioned in [14] but calculated at the limit of an infinite thickness: with in which  is the polar radius, while the parameters ] 1 , ] Expression ( 21) only depends on  44 and  66 via ] 3 .

Application to GaN
Gallium nitride is chosen to perform numerical applications because of its technological interest for electronic devices that can work under high power density and at relatively high temperature [4,5].The growth of GaN can be stabilized by epitaxial growth on bulky GaN substrates or a foreign substrate, for example, (111)Si, SiC, or sapphire [2,4,5,15].An investigation on GaN thin films is of particular interest because of its TDs that have anisotropic electrical properties regarding cathodoluminescence and photoluminescence.These properties are sensitive to the Burgers vectors and also to the TD directions [5].
The thickness of the GaN layer is assumed to be large enough to neglect the elastic interaction of the TD with the GaN/substrate interface.Different anisotropic elastic constants  of GaN are reported in several experimental or theoretical works.These constants are taken from [16][17][18][19][28][29][30] and are indicated in Table 1.GaN has the structure of wurtzite (space group P6 3 mc) with lattice parameters  = 0.319 nm and  = 0.519 nm [15].
Since the isotropic approximation is sometimes used in literature, column 7 indicates the calculated Young modulus  and the Poisson ratio ] relative to each  set. and ] are found from an averaging of the elastic properties of GaN according to the method proposed in [31].In this reference, a way is proposed to define an anisotropy factor coefficient  for any crystal symmetry.In the present work, using the  constants of [16][17][18][19][28][29][30], the  values are found between 1.27 and 1.78.These two values correspond to the  constants mentioned in [19] and [17], respectively.Previous transmission electron microscopy investigations, for example, [2,3,5,15], show that a large fraction of TDs cross the free surface with directions  parallel or not far from the c axis.They can be of edge-type, screw-type, or mixedtype.The local elastic equilibrium nearby the emerging point of a (a + c)-type TD can be treated from the image force theorem due to Lothe et al. [10].For instance, taking //a and a dislocation direction  in the prismatic plane (a, c), a systematic equilibrium angle  eq can be obtained.In Table 1, this angle is mentioned for each set of  constants given in [16][17][18][19][28][29][30] including isotropic approximation.The comparison of columns 8 and 9 shows a small anisotropic effect since  eq changes from averages equal to 2.5 ∘ (free surface, column 8) to 3.7 ∘ (no free surface, column 9).On the contrary, with the isotropic approximation, columns 10 and 11 indicate a strong effect of the free surface since  eq increases from averages equal to 4.7 ∘ (free surface) to 11.6 ∘ (no free surface).These last results demonstrate that the isotropic approximation is not reliable for GaN.
Other TDs in GaN are sometimes directed along the c axis, while other TDs have lines inclined noticeably with respect to c.For the screw case ( = 0) and b = c, the dislocation density b() is given by the expression (14).Taking the  constants indicated in [16], Figure 2 illustrates the sinusoidal shape of the two nonzero components of the dislocation density b(), that is,  1 () and  2 ().Our experience is that, for  = 0, the stress potential approach is more efficient to calculate the stress relaxation field than the integral formalism or the Fourier series.For  ̸ = 0, the Fourier approach and the integral formalism have similar performances.The numerical applications below, which use the elastic constants mentioned in line 2 Table 1, correspond to the geometry depicted in Figure 3.The TD is an inclined (a + c) screw dislocation placed in the  = 0 plane.Figure 4(a) describes the changes of the dislocation density in the plane  = 0.The component 3 is nonzero for this screw dislocation, Figure 4, conversely to the orientation  = 0.

Advances in Condensed Matter Physics
Figure 4(b) exhibits, as an example, the evolution of the  23 stress component along a straight line parallel to  and cutting the surface  = 0 at the point (, ) = (0, 2 nm).This inclined TD line is chosen because it gives some ideas about the results obtained from different approaches.All obtained curves converge far from the free surface but give different values close to the emerging point of the TD.The lower curve denoted "Aniso-∞" seems a bold line but is in fact the quasisuperposition of two curves, one calculated from the integral formalism and the other calculated from the Fourier series approach.For this latter approach, the calculation parameters were the following: period  = 30 nm, core size  = |b|.The squared geometry of the core has therefore no practical impact on the values of the stress at points placed about 2 or 3 nm from the dislocation axis.The number of harmonic terms (, ) chosen for summing the Fourier series was tested from visual convergence with the curve obtained from the integral approach.For points at 2 or 3 nm from the core, the number of harmonic terms should include a large number of terms for a good convergence ( and  can reach about a few hundreds).The same kind of difficulty is encountered with the integral formalism when it is applied to points close to the core.In the calculations, we adopted the integration domain ( − ) with  = 0.0001 rad.Table 1 shows that, with the other sets of  constants, the amplitudes of these curves change hardly.
The dashed curve denoted "Iso semi-∞" is calculated from the Yoffe-Shaibani-Hazzledine solution [6,7] and the upper curve denoted "Aniso-∞" assumes an infinite anisotropic crystal.The part of this latter curve, towards the negative , is not represented but is symmetrical with respect to  3 = 0.The bold line denoted "Aniso-∞" shows that both approaches presented in Sections 2 and 3 give the same result for the  23 stress and includes in particular the slight negative values nearby  3 = 0, that is, the region close to the dislocation core (a few nanometers).In this region, the stress is slightly negative, which is not the case in the isotropic approximation.This negative value nearby  3 = 0 is clearly due to the proximity to the TD core and not a calculation artefact generated by a lack of precision.Anisotropy effect is therefore far to be negligible nearby the free surface and the dislocation core.

Conclusions
In this work, the elastic field of an inclined TD has been investigated using different approaches.The first one, based on Fourier series and the concept of extended dislocation core, is developed in the present work from an appropriate combination of partial results obtained in [20,22].It proves to be a valuable alternative method to a well-known approach, the so-called integral formalism [10].The computation times required by both approaches are found roughly similar and depend on the wanted precision.There is a strong interest in using these two independent approaches to test the ability of elasticity to give the same values around the theoretical core.Both approaches converge beyond one nanometer from the theoretical core, as verified for some TDs in GaN.A third approach using stress potentials and assuming  = 0 is also numerically tested for the first time (at the author's knowledge).
When the anisotropic constants tend to those of an isotropic crystal, we verified numerically that the three above approaches converge to the same stress field, namely, that obtained from the Yoffe-Shaibani-Hazzledine expressions [6,7].Correction of misprints should be introduced in [7]; see [32].
For a dislocation normal to the basal free surface of a hexagonal crystal ( = 0), some general expressions of the integral formalism are simplified.They were tested positively in regions far from the free surface and close to the emerging core region of a TD.The numerical results exhibited in the four last columns of Table 1 and the curves in Figure 4(b) show clearly that the neglect of both elastic anisotropy and the presence of a free surface can be a confusion source in the interpretation of the angular equilibrium of a TD.
Of course, it is expected that the atomic structure of the dislocation core could influence the elastic field at distance from the core lower than a few Burgers vector moduli.According to our experience the assumption of a square section of the dislocation core has no sensible influence on the stress field beyond a few nanometers from the dislocation line provided that the core size  is less than the modulus of the Burgers vector.

A. Misprints in [26]
The following misprints can be noticed in Appendix C of [26].(iii) In expression  12 , the second denominator should be squared.

B. Misprint in [10]
In [10], the first expression in ( 14) should be corrected.All the trigonometric terms should be squared.

Figure 2 :
Figure 2: Case of the c screw TD oriented along the normal .Evolution of the dislocation density b( 1 ,  2 ,  3 ) in the basal plane  = 0 of GaN.

Table 1 :
Elastic constants of GaN.In anisotropic elasticity, the equilibrium angle  eq of the inclined (a + c) TD is small, weakly depending on the existence of a free surface (columns 8 and 9).It is not the case in isotropic elasticity (columns 10 and 11).Symbols A, S and I means Anisotropy, Surface and Isotropy.