On the importance of dislocation flow in continuum plasticity models for semiconductor materials

Journal of Crystal Growth 532 (2020) 125414

Contents lists available at ScienceDirect

Journal of Crystal Growth journal homepage: www.elsevier.com/locate/jcrysgro

On the importance of dislocation flow in continuum plasticity models for semiconductor materials

T

Binh Duong Nguyena, , Alexander M. Rauschb, Johannes Steinerc, Peter Wellmannc, Stefan Sandfelda ⁎

a

Chair of Micromechanical Materials Modelling, Institute of Mechanics and Fluid Dynamics, Technische Universität Bergakademie Freiberg (TUBAF), Lampadiusstraße 4, 09599 Freiberg, Germany b Chair of Materials Science and Engineering for Metals (WTM), Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU), Martensstraße 5, 91058 Erlangen, Germany c Crystal Growth Lab, Materials Department 6 (i-meet), Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU), Martensstraße 7, 91058 Erlangen, Germany

ARTICLE INFO

ABSTRACT

Communicated by Christopher Gourlay

Crystalline defects, such as dislocations, may have a strong negative influence on the performance and quality of electronic devices. Modelling the defect evolution of such systems helps to understand where defects occur, how they evolve and interact, Modelling also might allow to avoid – or at least to reduce – the number of line defects in such materials. A widely used continuum model for predicting the evolution of dislocation density in semiconductors is the Alexander-Haasen (AH) model (Alexander and Haasen, Solid state physics,1969) which describes the evolution of dislocation density through a local dislocation multiplication law; the motion of dislocations, i.e. dislocation flow, is not considered. Within this work, the underlying assumptions for the validity of the AH model are studied and compared to a more complex continuum model of dislocation dynamics (CDD) which explicitly considers dislocation fluxes. We use a simplified 2D scenario of a 4H-SiC crystal at the end of the growth as a benchmark system. Our comparisons show that considering the motion of dislocations can have a significant influence on the evolution and final distribution of dislocation density, which may strongly limit the applicability of the AH model. Based on the CDD model, we then investigate the cooling process and study the effect of the cooling rate on the resulting dislocation microstructure.

Keywords: A1. Computer simulation A1. Line defects A1. Stresses A1. Alexander-Haasen model B2. Semiconducting materials

1. Introduction Crystalline defects and imperfections may have a significant influence on the mechanical, chemical, or electronic properties of materials. This influence may be harmful or even useful, e.g., the dislocation, a linear defect, plays an essential role in the plastic deformations of crystalline materials – without which all crystalline materials would fail in a brittle manner. Volterra [2] introduced the dislocation as a purely mathematical concept, and in 1934, the existence of dislocations in the context of crystal plasticity was postulated by Orowan [3], Polanyi [4] and Taylor [5]. Subsequently, significant effort has been directed during the past decades towards experimental investigations as well as simulation methods for modelling and predicting dislocation structures. Lindroos [6] presented a new type of dislocation by using transmission electron microscopy technique to determine the Burgers vector of dislocations in Al-Mg alloys, Jenkinson and Lang [7] use X-ray diffraction topography showing dislocations in a single crystal of silicon, Kysar

⁎

et al. [8] investigate the dislocation structure formed in the nickel crystal after performing the wedge indentation and many other investigations as well. For simulation methods, this includes discrete dislocation dynamics [9–11], statistical mechanics of dislocations [12,13], strain gradient plasticity [14–16], continuum dislocation theory [17–20], continuum dislocation dynamics [21–25]. When it comes to the performance, reliability and behavior of many electronic devices the general design recommendations is to avoid dislocations since they might strongly deteriorate the functionality. Furthermore, the main focus of the above simulation approaches is on “after-growth” metals and alloys at room temperature. However, in the context of semiconductors it is crucial to understand the dislocation behavior also, e.g., during growth of crystals, at elevated temperatures, or during the cooling down process. Billig [26] and Tsivinsky [27] discovered that the dislocation density increases with the imposed temperature gradient and concluded that high thermal stress led to dislocation generation in the grown crystal. From this methodology,

Corresponding author. E-mail address: [email protected] (B.D. Nguyen).

https://doi.org/10.1016/j.jcrysgro.2019.125414 Received 5 October 2019; Received in revised form 27 November 2019; Accepted 6 December 2019 Available online 14 December 2019 0022-0248/ © 2019 Elsevier B.V. All rights reserved.

Journal of Crystal Growth 532 (2020) 125414

B.D. Nguyen, et al.

Jordan et al. [28] and later Miyazaki et al. [29] developed a simple phenomenological approach. They assume that the dislocation density is proportional to the total excess shear stress which is a total of the difference between the resolved shear stress and a critical resolved shear stress in each slip system. Later, Duseaux [30], Kobayashi and Iwaki [31] and Motakef [32] had used this approach to predict dislocation distribution in crystals grown from the melt. This approach was used to compute dislocation densities on a macroscopic level where the motion of dislocations is not visible. It could not verify a quantitative relation between the stress field and the dislocation density. The result can only be a very rough estimate of where and in which slip system dislocations may be generated during crystal growth. This already showed at an early point of time a strong need for a more advanced model that ideally include both dislocation motion as well as multiplication. Alexander and Haasen [1] derived a single slip model, the so-called Alexander-Haasen (AH) model for describing the evolution of the mobile dislocation density and plastic deformation on one primary slip system. Yonenaga and Sumino [33], Suezawa et al. [34] developed a one-dimensional material response function for silicon, based on Haasen’s work. It is an elastic-viscoplastic type of constitutive equation which includes dislocation densities as internal state variables. It represents a significant advancement to predict the quantity of dislocation density as a result of temperature and strain rate. Base on that, Tsai [35] has developed a nonlinear finite element formulation for the prediction of dislocation densities and residual stresses during crystal growth which was applied for studying a variety of crystals such as GaAs, InP, Si, SiC [e.g. 36–38]. Recent work focusing on the dislocation evolution during crystal growth was done by Gao et al. [39], Ide et al. [40]. They were the first to compare dislocation density distributions from the AH model with experimental data for Si to obtain a more indepth insight into the validity of the model. However, the propagation of dislocations is still not considered in these rather phenomenological approaches which might limit their predictive power. In this work, we take a “bottom-top” approach towards plasticity which builds upon the motion and interaction of discrete dislocations, represented by continuous fields in a non-local framework of dislocation fluxes. There, the motion of dislocations gives rise to dislocation structures that can not be described by local phenomenological models. Our model is the combination of AH model and the Groma model [41], which does not differentiate between mobile and immobile dislocations but instead sorts them into edge dislocations of positive and negative orientation. This model is systematically derived from coarse-graining discrete dislocations and their collective dynamics. While this model is still not as accurate or detailed as a discrete dislocation dynamics model, it has the efficiency of a continuum representation where even large numbers of dislocations are just represented through a number: the density. In the following, we start by introducing the theoretical formulation of the AH model, followed by the flux-based Groma model. Then, simulation results are presented comparing the AH model with the Groma model, which we extended for comparison purposes such that it accounts for dislocation multiplication in the same manner as the AH model. Furthermore, we investigate the influence of different cooling rates. We conclude by a discussion and give recommendations when dislocation fluxes should be explicitly considered in a model and when a less accurate model might be sufficient.

in a slip system with the slip plane normal vector n and the slip direction s

y

res

= ni

(1)

ij sj

= aGb

(2)

m

with a constant a, the shear modulus G, the magnitude of the burgers vector b and the mobile dislocation density m . This leads to an effective stress acting on dislocations

= |

eff

res |

with

y

· =

· for · > 0 . 0 for · 0

(3)

The average velocity of a dislocation is based on a dislocation core diffusion model [44]: m

eff

v = v0

Q sign( kB T

exp

0

res )

(4)

with a constant v0 , the stress exponential factor m, the activation energy Q, the Boltzmann constant kB , and the absolute temperature T. The evolution of m is described by a local dislocation multiplication law that originally was proposed by Johnston and Gilman [45] and refined by Peissker et al. [46],

=K

m

eff 0

m |v|

(5)

where K and are two constants. Finally, the Orowan equation [47] describes the evolution of the plastic shear strain pl in the slip system due to dislocation motion pl

=

m bv.

(6)

2.2. Groma model including nucleation Contrary to the AH model, the Groma model [41] is a non-local model that describes the flux of dislocation density by a set of transport equations. As for the AH model we also use the Groma model only in a single slip setting. Dislocation densities in the Groma model consist of straight edge dislocations which may have two orientations (a negative edge dislocation has opposite line direction as that of a positive edge). They pierce a two-dimensional plane which is why they can essentially be represented as point-like “particles”. Then, Burgers vector and slip plane normal are aligned with the x- and z-axis of the coordinate system (see Fig. 1a). Accordingly, the positive and negative dislocation densities, + and , are defined as the number of positive or negative dislocations per averaging volume. Then evolution equations are transport equations formulated for the total and signed geometrically necessary dislocation densities, tot = + + and gnd = + , as tot gnd

= =

x ( gnd v ) x ( tot v ).

+S

(7) (8)

where S is a source term governing dislocation nucleation. This is a phenomenological representation of the Frank-Read source mechanism, which is responsible for dislocations being generated on the chosen slip plane. As stated in [48] for a 2D system of discrete edge dislocations, when the applied shear stress at the source exceeds the critical value during a critical time duration, the newly nucleated configuration is formed by bowing out of the dislocation segment until a full loop is formed. This process repeats as long as the shear stress exists larger or equal the critical value. Thus, dislocations are being generated and multiplied. The continuous source term is a spatio-temporal average of this process. The velocity law and the formulation of stresses are chosen

In the following, we introduce the theoretical background of the AH model [1] as used in our simulation. We then briefly introduce the Groma model [41] and explain our approach for considering dislocation nucleation. 2.1. Alexander-Haasen model res

or

where is the stress tensor in the local slip system coordinate system which has the components ij . Due to friction-type dislocation interactions res is offset by a yield stress [42,43]

2. Theory

Dislocation motion is driven by the resolved shear stress

= n · ·s

res

acting 2

Journal of Crystal Growth 532 (2020) 125414

B.D. Nguyen, et al.

Fig. 1. (a) Edge dislocation with negative (left) and positive (right) line direction in 2D space with a slip direction and slip normal vector aligned with the x- and zdirection. (b) Cut through a cylindrical crystal (indicated by the blue rectangle) with height h and width 2r representing the considered 2D crystal geometry, (c) the temperature field, and (d) the resulting initial resolved shear stress res in the basal plane for this temperature field. The red dashed line in (b) indicates the data position used for line plots. At the points C1 and C2 mechanical boundary conditions are applied as explained in the text.

as in the AH model. According to [41] the source term is proportional to the product of the stress (here the effective stress eff ) and the evolution of the plastic strain pl ,

S=C

eff tot v

boundary conditions, the top of the crystal is fixed in the vertical direction and can freely move in the horizontal direction: for the points shown in Fig. 1b it is at C1: u x = uz = 0 and at C2: uz = 0 with u x and uz being the displacements in x- and z-direction. Thermal insulation boundary conditions are applied at the crystal boundaries. An axialsymmetrical temperature field which does not change over time will be prescribed (see Fig. 1c) during the simulations. The bottom of the crystal has a typical process temperature of roughly 2300 °C . The temperature gradient between the bottom and the top is about 100 °C in this model system and causes thermal stress inside the crystal. The resulting resolved shear stresses are shown in Fig. 1d and are in the range of 2 MPa . Compared with simulations results from literature resolved shear stresses in the range of 2 MPa for a PVT process are reasonable for hexagonal SiC [50]. Note, that in reality the stresses also depend on details of the interface between the crystal seed and the holder [51]. The critical resolved shear stress in the basal plane for observing dislocation activity in 4H-SiC in compression test is below 2 MPa if the data from [52–54] are extrapolated to a temperature of 2500 K. Thus, these stresses are comparably low as compared to the thermal stresses. Therefore, neither the AH nor the Groma model contain any critical stress that would be required to activate dislocation multiplication or dislocation propagation.

(9)

with a constant C. This equation has some similarity with (5) of the AH model. In the following we only want to study the difference due to the transport character of the evolution equations. Therefore, we use the source term from the AH model as well in the Groma model by assuming that tot m . The subscript “m ” in the AH model indicates that only mobile and no sessile dislocations are considered; the Groma model does not have to differentiated between the two “types” of dislocations. However, since in our model no, e.g. junction formation can take place, all dislocations could be considered as mobile. The evolution of the plastic shear strain is calculated based on the Orowan Eq. (6) (replacing m by tot ). 3. Results 3.1. Simulation setup and geometry As an example, hexagonal SiC (4H) was considered that was grown by a typical physical vapor transport (PVT) process. For simplicity, we assume an idealized two-dimensional SiC crystal as shown in Fig. 1b with a typical temperature distribution, but we do not consider the crystal growth process itself. In order to relate modelling to the experimental data from [49] the dimensions of PVT growth SiC crystal are set to 45 mm and 30 mm for the height and radius respectively. The 2D system can be thought of as a vertical slice through a cylindrical SiC crystal. The dashed line at z = 9 mm in Fig. 1b defines the line used for some plots later on. The basal plane is considered as the primary and only active slip system. The plane normal points along the z-direction such that slip happens parallel to the x-direction. As mechanical

3.2. Material parameter The parameters for the AH model are taken from [55] and are summarized in the appendix in Table A.1. The temperature dependent shear modulus is taken as the elastic constant C44 , see Table A.2. Instead of using two stress terms (long-range internal stress and short-range stress) as in [55] for the computing effective stress, we only consider y as in Eq. (3). Accordingly, the value for the factor a is a sum of the factors in both stress terms by assuming only one active slip system. The temperature dependent elastic constants for hexagonal SiC are 3

Journal of Crystal Growth 532 (2020) 125414

B.D. Nguyen, et al.

extrapolated from the temperature dependent values by Li and Bradt [56] (see Table A.2). According to Zhang et al. [50] an isotropic temperature dependent thermal expansion coefficient (T ) is considered:

(T ) =

5.54·10 2.78·10

4

6

K 1exp

K

1exp

(

(

T 3358.8 K T 42.3 K

4.1.1. Analysis of the Alexander-Haasen model To understand the evolution of the mobile dislocation density m for the AH model, Fig. 2e shows the resolved shear stress res field for the basal slip system at t = 3000 s and the evolution of the absolute value of the resolved stress | res |. Initially, m will increase at the four stress extrema (compare Fig. 2c and for a later stage also Fig. 2a). Due to the local increase of the dislocation density, the stresses in these regions decrease as a consequence of plastic relaxation (compare Fig. 2e). This leads to a additional redistribution of stresses as can be observed in the shape change (Fig. 1d vs Fig. 2e). In particular, during the evolution the stress increases sharply at the borders of the plastic relaxation zone (compare Fig. 2g). Accordingly, the dislocation density increases along a sharp line (see Fig. 2a). As a result of this stress transfer, a substantial increase in dislocation density by nucleation of new dislocations takes place in regions where originally the initial stress was low. The evolution of the plastic shear strain follows that of the dislocation density (Fig. 2i and Fig. 2k), as would be expected (combining (5) and (6) gives pl m b ). Therefore, the shapes of the non-zero regions are very similar as compared to the non-zero density patches. For the AH model the plastic strain only increases along with the multiplication of dislocation density: “high density = high plastic strain”.

)

) + 8.83·10

6

K 1.

(10)

3.3. Initial values and boundary conditions for the dislocation problem The initial dislocation densities ( m or tot ) need to be set to a non zero value, because the source term in the AH model accounts only for dislocation multiplication and not for nucleation of new dislocation sources. Accordingly, if no dislocation density is present at the beginning, the dislocation density will not change. Thus, the initial mobile dislocation density in the AH model and the total dislocation density in the Groma were set to a small value of 1 cm 2 . Additionally, the initial GND density was set to zero for the Groma model. In the Groma model, a mathematical boundary value problem must be solved, i.e., we have to prescribe what happens when dislocations reach the surface. Here, we assume impenetrable boundaries where dislocations can not leave the crystal. In physical terms, this behavior could be the result of an oxide layer/different polycrystalline structure.

4.1.2. Analysis of the Groma model combined with the AH source term Combining the formulation of source term for nucleation from the AH model with the pure transport equations, as in (7) and (8), shows a density and plastic strain evolution where only few similarities with the AH model can be found. In Fig. 2b one can observe that right from the very beginning the distinct density accumulation at the locations of the stress extrema is absent (compare positions marked by label A and B in Fig. 2c and Fig. 2d). At the same time dislocation accumulations occur close to the surface and close to the crystal axis. The reason for this behavior is that dislocations are nucleated where the stress is highest, but then move away, driven by shear stresses. They accumulated where the shear stress is low (see Fig. 2f and Fig. 2h) and form quasi-stationary pile-up structures. This is the case both for the region near the axis (label C in Fig. 2d) as well as for the surface (label D) and ultimately leads to a depleted zone where in the AH model the density is highest. Fig. 3 shows the evolution of dislocation density over time. The dislocation density propagates away from the stress hot spots (compare Fig. 1d) along the basal plane in positive or negative x-direction towards either the surface or the middle axis (see Fig. 2b and Fig. 2d); surfaces are impenetrable for dislocations. Consequently, dislocation density will pile up there. According to the Groma model only positive and negative edge dislocation are present. At the beginning the density for positive and negative edge dislocations is the same everywhere = 0 ). Depending on the sign for the stress hot spots one ( gnd = + kind of either positive or negative edge dislocation will travel in positive or negative x-direction. For example, the negative hot spot on the bottom left is considered (see Fig. 2f). Here, negative dislocations will move towards the middle axis, positive will move towards the surface. This can be seen in Fig. 4, where the geometrically necessary gnd is plotted. A positive value represents an excess of positive and a negative value an excess of negative dislocations. Thus, if dislocations moving towards the middle axis would meet from opposite stress hot spots, they could not annihilate due to the same sign. Fig. 4 shows the geometrically necessary dislocation density gnd for the Groma case with multiplication. Here it can be seen that the dislocation density towards the middle axis consists of an excess of the opposite kind of dislocations that are moving to the corresponding surface. Deviating from the AH model, the plastic shear strain in the Groma model is caused by dislocation motion. In Fig. 2j and Fig. 2l we can see that the plastic shear strain is high in the regions which were swept by dislocation.

3.4. Numerical details For solving the 2D problem the commercial finite element software COMSOL Multiphysics® is used. Gallien et al. [57] showed the compatibility of solving the AH equations with COMSOL Multiphysics®. For the transport equations in the Groma model, a small artificial diffusion term is added to minimize oscillations. The finite element mesh is the same for all simulations. The mesh was chosen fine enough, so that further refinement does not change the results. Quadratic shape functions are used in all cases. The backward differentiation formula (BDF) method is chosen for the time integration and the smallest time step in this case is 10 s. 4. Results and discussion Two different simulation setups are discussed subsequently: we start by studying the occurring dislocation structure during a relaxation process at constant temperature and will then investigate the influence of an additional cooling-down process. 4.1. Evolution of dislocation densities at constant temperature The following comparisons consider the AH model (pure multiplication model) and compare it to the Groma model. To make the models comparable, the Groma model also includes the same source term for dislocations multiplication as the AH model. Simulation results of the evolution of dislocation density at an early stage along the red, dashed line defined in Fig. 2a and Fig. 2b are shown in Fig. 2c and Fig. 2d. The corresponding 2D density distributions for the systems at a later time step (t = 3000 s ) are shown in Fig. 2a and Fig. 2b. We observe a larger difference – not only in terms of magnitudes of the density values but also in the qualitative shape of the plots: The AH model shows a growing region of increasing density, while the Groma model exhibits two regions where the density has maxima (one at the surface, x 30 mm and one closer to the crystal axis, x = 0 mm). These differences even become more extreme when the simulation progresses. In order to understand which details of the underlying physics are represented by which model we will now analyze both simulations in more detail.

4.1.3. Discussion Comparing the resulting dislocation densities for the local AH model with the nonlocal formulation (AH source term together with the 4

Journal of Crystal Growth 532 (2020) 125414

B.D. Nguyen, et al.

Fig. 2. Dislocation density m (AH) or tot (Groma) in (a) and (b), resolved shear stress res in (e) and (f) and plastic shear strain in (i) and (j) after 3000 s for the cases AH model and Groma model with multiplication; and plot along the cut line for several times (red dashed line) respectively in (c) and (d), (g) and (h), (k) and (l). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Groma model) shows that dislocation propagation has a non negligible influence. The need for an explicit propagation law can also be seen based on the following estimate: For the pure AH setup, the velocity field is obtainable for every time step (here 1 s ). Thus, by choosing an arbitrary starting point for an edge dislocation the distance a dislocation travels along the slip plane for every time step can be calculated by integrating the average velocity at the current position for every step.

Note, that in doing so one implicitly assumes that dislocations are not immobilized. The results for several times for a positive and negative dislocation with a starting point lying in the middle of the stress hot spot on the bottom left position of the system can be seen in Fig. 4. This estimate also predicts that dislocations will propagate to the surface and towards the middle axis. If the edge dislocation will reach the surface it would create a surface step and “leave” the material. As for the 5

Journal of Crystal Growth 532 (2020) 125414

B.D. Nguyen, et al.

Fig. 3. Evolution of dislocation density with the increasing time step in Groma model.

simulations this was neglected here. Basically this gives a good approximation for the results obtained for the Groma case with multiplication (compare Fig. 4). However, using the source term from the AH model in the Groma model is a rough approximation by itself. All the parameters for the AH model are determined from a set of stress-strains curves at different temperatures [58]. It is not clear, whether this parameters are suitable for arbitrary deformation conditions. Furthermore, in terms of imperfections in the defect geometry and crystallography, only edge dislocations in the basal plane - so called basal plane dislocations(BPDs) are considered in the simulations whereas the AH model does not differ between different dislocation types in real system such as threading edge dislocations(TEDs), threading screw dislocations(TSDs), and superscrew dislocations - called micropipes (MPs), and the interaction between them which may lead to the formation of low angle grain boundaries(LAGBs). Thus, for more reliable results, a more consistent model would be required. Nevertheless, Ha et al. [51] mentioned the possibility that dislocations may propagate to the center of the crystal during crystal growth, pile up and causing bending of the crystal. The accumulation of dislocations towards the center and the periphery of PVT grown 6H-SiC crystals can be found in the basal plane for etched wafers from Sakwe et al. [49] (see Fig. 5c) as the results in Fig. 2b also

indicate for the Groma model with multiplication. To obtain an impression how the results for the dislocation densities would look like in the xy-plane (basal plane) instead of the considered xz-plane in the simulations, the 2D dislocation densities for the AH and Groma with multiplication case are first integrated along the y-direction of the crystal for the cut line defined in Fig. 1. This results in a projection of the 2D system onto the x-axis. The 1D data can now be revolved around the central axis of a virtual cylindrical crystal. The resulting output mimics the dislocation distribution in an on-axis SiC wafer (see Fig. 5). Note, that there are three symmetry equivalent slip directions in the basal plane, whereby here the 1D data from only one slip direction continuously extended. Also one example from [49] showing qualitatively the dislocation density distribution (darker regions correspond to higher dislocation densities) from an SiC wafer can be seen in Fig. 5. Here, the density is highest in a large region around the centre and at the periphery of the wafer. The hexagonal symmetry in the basal plane is observable due to darker spots in the brighter ring between the high dislocation densities in the middle and the periphery. In the AH case, only the dislocation density is high everywhere except for the centre. As mentioned before, the results from the Groma system indicate that the dislocation density at the periphery of the wafer shown in the experimental data may be due to dislocations moving to the

Fig. 4. Estimation for the travel distance for a positive and negative dislocation in the basal plane originated at the location of maximum resolved shear stress at the beginning for several times. The 2D data shows the patches of GND dislocation, respectively. 6

Journal of Crystal Growth 532 (2020) 125414

B.D. Nguyen, et al.

Fig. 5. Revolution plots for AH, Groma with multiplication and an etched wafer showing qualitatively the dislocation distribution in the basal plane of a nominally undoped n-type 6H-SiC crystal [49] (darker regions correspond to higher dislocation density).

surface of the crystal. Dislocation density towards the middle axis can also be found for the Groma case. For arbitrary materials and growth or cooling processes, it is not possible to say whether dislocation transport has a significant influence on the evolution of dislocation density in general. The dominance of either multiplication or propagation strongly depends on the ratio of how fast dislocation density increases locally compared to how fast they can move away from their origin. A fast local increase leads to plastic relaxation. Also, the higher the dislocation density, the lower will be the mean free path of dislocation, because of faster relaxation and increase of the yield stress y . Thus, dislocation transport may not be relevant on a macroscopic scale. However, if the increase in dislocation density is slow, dislocations may have enough time to propagate over distances in the range of the crystal dimensions. Also Gao et al. mentioned in [39] that for Si the explicit consideration of propagation might be important: They qualitatively compared simulation results from an extended AH model with experimental data. The degree of accordance of the simulations with the experiments seemed to be dependent on the crystal orientation during crystal growth. Their conclusion was that the discrepancy was due to the lack of an explicit dislocation propagation law which is in line with our findings.

that in real system the temperature gradient distribution is different, that it is changed at each time step and the material is not perfectly homogeneous. Fig. 7 shows snapshots of the resolved shear stress distribution. Due to the large difference between the stress at 2500 K and stress at room temperature, two different ranges of the colour bar are used here for more convenient observation. From t = 0 s to t = 500 s, because of the constant temperature, the maximum of stress distribution is around 2 MPa. Immediately after the end of the first regime in case 1 the stress suddenly drop in response to the rapid cooling-down and reaches a value of around 0.1 MPa at t = 510 s, while stress hot spot in case 2 and case 3 still at high values. At t = 2000 s the temperature in case 2 now reaches room temperature, and the stress also drops to around 0.1 MPa while in case 3 stress still remains high value. At t = 4000 s, we can see that the resolved shear stress in all 3 cases is now in the same range with the maximum value of around 0.1 MPa. Although the thermal stresses are in the same range in three cases, the stress from dislocations of negative sign is higher in case 2 and highest in case 3 as compared to case 1. During relaxation, the dislocations propagate away from the stress hot spot are shown in Fig. 3. A very high stress concentration causes a high nucleation rate, therefore, in case 1 with rapid cooling, dislocation density remain constant after t = 500 s, while dislocations still multiply until t = 2000 s in case 2 and t = 3000 s in case 3 as can be seen in Fig. 6b. After that, at room temperature, dislocations stop nucleating and the density remains at constant value. Due to plastic relaxation in the models (see Eq. 6) during the simulations, stresses are relieved by dislocation rearrangement until eventually a steady state (at t = 4000 s) is reached. Dislocation density distributions at this state are shown in Fig. 8a together with the GND density distribution in Fig. 8b. The crystal which is cooled down over a longer duration contains a much higher density of dislocation and the crystal which is rapidly cooled down results in the lowest dislocation density in comparison to the other crystals (Fig. 6b and Fig. 8a). These results explain the experimental results [59] and show good qualitative agreement with.

4.2. Investigation of a cooling process based on the Groma model While the above section investigated the dislocation evolution during relaxation of the crystal we will now combine two different stages: relaxation from t = 0 s to t = 500 s, followed by a cooling phase of up to 4000 s where the temperature decreased down to room temperature ( 20 °C). For the cooling phase three different temperature rates are studies (compare Fig. 6a): an instantaneous temperature drop (case 1) and a decrease to room temperature by a rate of approximately 1.47 K/s (case 2) and 0.88 K/s (case 3). Subsequently, we will use the Groma model to investigate how dislocation multiplication and propagation is affected by the cooling. In this first cooling process approach, for simplification, we assume that the temperature gradient distribution remains constant during cooling down although the fact

Fig. 6. Cool down temperature (a) and integration of total dislocation density (b) respectively. 7

Journal of Crystal Growth 532 (2020) 125414

B.D. Nguyen, et al.

Fig. 7. Resolved shear stress from t = 0 s to t = 4000 s for the right bottom quadrant the two-dimensional cross-section shown in Fig. 1b. Two different color scales are used to visualize the largely different range of values.

5. Conclusion

dislocation model would give a more in-depth insight into how crucial the neglection of an explicit dislocation propagation law in the AH model is.

In summary, the comparison of the AH model with the Groma model with dislocation multiplication showed that dislocation transport may strongly influence the evolution and final distribution of the dislocation density. The degree of this influence generally depends on the rate of the local increase of dislocation density and how fast dislocation density can propagate away from their origin. We also show that the rate of the cooling process has a significant effect on dislocation multiplication. However, the simplified 2D SiC example is still far away from realistic crystal growth processes. The time in the simulations presented here for several minutes is only a small excerpt of a real growth or cooling process that would last for hours or days. The longer the cooling process is, the larger the dislocation density is in the crystal. A 3D crystal growth framework with a more general continuum

Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgement Financial support of the Deutsche Forschungsgemeinschaft (DFG) under the contract numbers WE2107-15 and SA2292-6 is greatly acknowledged.

Fig. 8. Total dislocation density (a) and GND density (b) distribution at t = 4000 s. 8

Journal of Crystal Growth 532 (2020) 125414

B.D. Nguyen, et al.

Appendix A. Parameters for the simulations See Table A.1 and A.2.

Table A.1 Parameters for the AH model for SiC. Parameter

Value −10

m 3.073·10 8.5·10−15 ms−1 2.0·10−5 m−1 3.3 eV for T > 1000 °C 8.617·10−5 eV K−1 1.1 for T > 1000 °C 2.8 0.425

Magnitude of Burgers vector b (< 1 1 2¯ 0>) Velocity prefactor v0 Multiplication factor K Activation energy Q Boltzmann constant kB Stress exponent Stress exponent m Prefactor a

Table A.2 Anisotropic temperature dependent elastic constants for hexagonal SiC (Voigt notation). Parameter

Value

C11

0.025 GPa K 1· T + 486.6 GPa K

C12

0.011 GPa K 1· T + 101.3 GPa K

C13

0.011 GPa K 1· T + 59.02 GPa K

C33

0.025 GPa K 1· T + 528.9 GPa K

C44

0.007 GPa K 1· T + 150.3 GPa K

References

dislocations, Int. J. Plast. 12 (5) (1996) 611–627. [20] K.C. Le, Thermodynamic dislocation theory for non-uniform plastic deformations, J. Mech. Phys. Solids 111 (2018) 157–169. [21] I. Groma, F.F. Csikor, M. Zaiser, Spatial correlations and higher-order gradient terms in a continuum description of dislocation dynamics, Acta Mater. 51 (5) (2003) 1271–1281. [22] T. Hochrainer, S. Sandfeld, M. Zaiser, P. Gumbsch, Continuum dislocation dynamics: towards a physical theory of crystal plasticity, J. Mech. Phys. Solids 63 (2014) 167–178. [23] S. Sandfeld, T. Hochrainer, P. Gumbsch, M. Zaiser, Numerical implementation of a 3d continuum theory of dislocation dynamics and application to micro-bending, Phil. Mag. 90 (27–28) (2010) 3697–3728. [24] S. Sandfeld, T. Hochrainer, M. Zaiser, P. Gumbsch, Continuum modeling of dislocation plasticity: theory, numerical implementation, and validation by discrete dislocation simulations, J. Mater. Res. 26 (5) (2011) 623–632. [25] S. Sandfeld, M. Monavari, M. Zaiser, From systems of discrete dislocations to a continuous field description: stresses and averaging aspects, Modell. Simul. Mater. Sci. Eng. 21 (8) (2013) 085006. [26] E. Billig, Some defects in crystals grown from the melt-i. defects caused by thermal stresses, Proc. R. Soc. Lond. Ser. A. Math. Phys. Sci. 235 (1200) (1956) 37–55. [27] S.V. Tsivinsky, Dislocation density in pure crystals grown from melts, Kristall und Technik 10 (1) (1975) 5–35. [28] A.S. Jordan, A.R. Von Neida, R. Caruso, The theory and practice of dislocation reduction in gaas and inp, J. Cryst. Growth 70 (1–2) (1984) 555–573. [29] N. Miyazaki, H. Uchida, S. Hagihara, T. Munakata, T. Fukuda, Thermal stress analysis of bulk single crystal during czochralski growth (comparison between anisotropic analysis and isotropic analysis), J. Cryst. Growth 113 (1–2) (1991) 227–241. [30] M. Duseaux, Temperature profile and thermal stress calculations in gaas crystals growing from the melt, J. Cryst. Growth 61 (3) (1983) 576–590. [31] N. Kobayashi, T. Iwaki, A thermoelastic analysis of the thermal stress produced in a semi-infinite cylindrical single crystal during the czochralski growth, J. Cryst. Growth 73 (1) (1985) 96–110. [32] S. Motakef, Fundamental considerations in creep-based determination of dislocation density in semiconductors grown from the melt, J. Cryst. Growth 114 (1–2) (1991) 47–58. [33] I. Yonenaga, K. Sumino, Dislocation dynamics in the plastic deformation of silicon crystals i. Experiments, Phys. Status Solidi (A) 50 (2) (1978) 685–693. [34] M. Suezawa, K. Sumino, I. Yonenaga, Dislocation dynamics in the plastic deformation of silicon crystals. II. Theoretical analysis of experimental results, Phys. Status Solidi (a) 51 (1) (1979) 217–226. [35] C.T. Tsai, On the finite element modeling of dislocation dynamics during semiconductor crystal growth, J. Cryst. Growth 113 (3–4) (1991) 499–507. [36] J. Völkl, G. Müller, A new model for the calculation of dislocation formation in semiconductor melt growth by taking into account the dynamics of plastic

[1] H. Alexander, P. Haasen, Dislocations and plastic flow in the diamond structure, Solid State Physics, vol. 22, Elsevier, 1969, pp. 27–158. [2] V. Volterra, Sur l’équilibre des corps élastiques multiplement connexes, Annales scientifiques de l’École normale supérieure 24 (1907) 401–517. [3] E. Orowan, Zur kristallplastizität. ii, Z. Phys. 89 (9–10) (1934) 614–633. [4] M. Polanyi, Über eine art gitterstörung, die einen kristall plastisch machen könnte, Z. Phys. 89 (9–10) (1934) 660–664. [5] G.I. Taylor, The mechanism of plastic deformation of crystals. Part I. Theoretical, Proc. R. Soc. Lond. Ser. A Contain. Pap. Math. Phys. Charact. 145 (855) (1934) 362–387. [6] V.K. Lindroos, Knitting of a dislocation network with selective interaction, Phil. Mag. 24 (189) (1971) 709–712. [7] A.E. Jenkinson, A.R. Lang, X-ray diffraction topographic studies of dislocations in floating-zone grown silicon, in: Conference on Direct Observation of Imperfections in Crystals, 1962, pp. 471–495. [8] J.W. Kysar, Y. Saito, M.S. Oztop, D. Lee, W.T. Huh, Experimental lower bounds on geometrically necessary dislocation density, Int. J. Plast. 26 (8) (2010) 1097–1123. [9] H.D. Espinosa, M. Panico, S. Berbenni, K.W. Schwarz, Discrete dislocation dynamics simulations to interpret plasticity size and surface effects in freestanding fcc thin films, Int. J. Plast. 22 (11) (2006) 2091–2117. [10] C. Motz, D. Weygand, J. Senger, P. Gumbsch, Initial dislocation structures in 3-d discrete dislocation dynamics and their influence on microscale plasticity, Acta Mater. 57 (6) (2009) 1744–1754. [11] D. Mordehai, E. Clouet, M. Fivel, M. Verdier, Introducing dislocation climb by bulk diffusion in discrete dislocation dynamics, Phil. Mag. 88 (6) (2008) 899–925. [12] S. Yefimov, I. Groma, E. Van der Giessen, A comparison of a statistical-mechanics based plasticity model with discrete dislocation plasticity calculations, J. Mech. Phys. Solids 52 (2) (2004) 279–300. [13] O. Kapetanou, V. Koutsos, E. Theotokoglou, D. Weygand, M. Zaiser, Statistical analysis and stochastic dislocation-based modeling of microplasticity, J. Mech. Behav. Mater. 24 (3–4) (2015) 105–113. [14] N.A. Fleck, G.M. Muller, M.F. Ashby, J.W. Hutchinson, Strain gradient plasticity: theory and experiment, Acta Metall. Mater. 42 (2) (1994) 475–487. [15] S. Wulfinghoff, T. Böhlke, Gradient plasticity for single crystals, PAMM 10 (1) (2010) 351–352. [16] L.P. Kubin, A. Mortensen, Geometrically necessary dislocations and strain-gradient plasticity: a few critical issues, Scr. Mater. 48 (2) (2003) 119–125. [17] K. Kondo, On the geometrical and physical foundations of the theory of yielding, Proc. 2nd Japan Nat. Congr. Applied Mechanics, vol. 2, 1952, pp. 41–47. [18] J.F. Nye, Some geometrical relations in dislocated crystals, Acta Metall. 1 (2) (1953) 153–162. [19] K.C. Le, H. Stumpf, A model of elastoplastic bodies with continuously distributed

9

Journal of Crystal Growth 532 (2020) 125414

B.D. Nguyen, et al. deformation, J. Cryst. Growth 97 (1) (1989) 136–145. [37] D. Maroudas, R.A. Brown, On the prediction of dislocation formation in semiconductor crystals grown from the melt: analysis of the haasen model for plastic deformation dynamics, J. Cryst. Growth 108 (1–2) (1991) 399–415. [38] B. Gao, K. Kakimoto, Numerical investigation of the influence of cooling flux on the generation of dislocations in cylindrical mono-like silicon growth, J. Cryst. Growth 384 (2013) 13–20. [39] B. Gao, K. Jiptner, S. Nakano, H. Harada, Y. Miyamura, T. Sekiguchi, K. Kakimoto, Applicability of the three-dimensional alexander-haasen model for the analysis of dislocation distributions in single-crystal silicon, J. Cryst. Growth 411 (2015) 49–55. [40] T. Ide, H. Harada, Y. Miyamura, M. Imai, S. Nakano, K. Kakimoto, Relationship between dislocation density and oxygen concentration in silicon crystals during directional solidification, Crystals 8 (6) (2018) 244. [41] I. Groma, P. Balogh, Investigation of dislocation pattern formation in a two-dimensional self-consistent field approximation, Acta Mater. 47 (13) (1999) 3647–3654. [42] G.I. Taylor, The mechanism of plastic deformation of crystals. Part I. Theoretical, Proc. R. Soc. Lond. Ser. A Contain. Pap. Math. Phys. Charact. 145 (855) (1934) 362–387. [43] G.I. Taylor, The mechanism of plastic deformation of crystals. Part II. Comparison with observations, Proc. R. Soc. Lond. Ser. A Contain. Pap. Math. Phys. Charact. 145 (855) (1934) 388–404. [44] P. Haasen, On the plasticity of germanium and indium antimonide, Acta Metall. 5 (10) (1957) 598–599. [45] W.G. Johnston, J.J. Gilman, Dislocation velocities, dislocation densities, and plastic flow in lithium fluoride crystals, J. Appl. Phys. 30 (2) (1959) 129–144. [46] E. Peissker, P. Haasen, H. Alexander, Anisotropic plastic deformation of indium antimonide, Phil. Mag. 7 (80) (1962) 1279–1303. [47] E. Orowan, Problems of plastic gliding, Proc. Phys. Soc. 52 (1) (1940) 8. [48] E. Van der Giessen, A. Needleman, Discrete dislocation plasticity: a simple planar model, Modell. Simul. Mater. Sci. Eng. 3 (5) (1995) 689.

[49] S.A. Sakwe, M. Stockmeier, P. Hens, R. Müller, D. Queren, U. Kunecke, K. Konias, R. Hock, A. Magerl, M. Pons, A. Winnacker, P. Wellmann, Bulk growth of SiC - review on advances of SiC vapor growth for improved doping and systematic study on dislocation evolution, Phys. Status Solidi (B) 245(7) (2008) 1239–1256. ISSN 15213951. [50] Z. Zhang, J. Lu, Q. Chen, V. Prasad, Thermoelastic stresses in sic single crystals grown by the physical vapor transport method, Acta. Mech. Sin. 22 (1) (2006) 40–45. [51] S. Ha, G.S. Rohrer, M. Skowronski, V.D. Heydemann, and D.W. Snyder. Plastic deformation and residual stresses in sic boules grown by pvt, in: Materials Science Forum, vol. 338, Transtec Publications, 1999, 2000, pp. 67–70. [52] A.V. Samant, P. Pirouz, Activation parameters for dislocation glide in α-SiC. Int. J. Refractory Met. Hard Mater. 16(4) (1998) 277–289. ISSN 0263-4368. [53] P. Pirouz, M. Zhang, J.L. Demenet, H.M. Hobgood, Yield and fracture properties of the wide band-gap semiconductor 4H-SiC, J. Appl. Phys. 93 (6) (2003) 3279–3290. [54] A. Lara, A. Muñoz, M. Castillo-Rodríguez, and M. Domínguez-Rodríguez. Plastic behaviour of 4H-SiC single crystals deformed at temperatures between 800 and 1300°C. Ceram. Int. 38(2) (2012) 1381–1390. ISSN 0272–8842. [55] B. Gao, K. Kakimoto, Three-dimensional modeling of basal plane dislocations in 4HSiC single crystals grown by the physical vapor transport method, Cryst. Growth Des. 14 (3) (2014) 1272–1278. [56] Z. Li, R.C. Bradt, The single crystal elastic constants of hexagonal SiC to 1000°C, Int. J. High Technol. Ceram. 4 (1) (1988) 1–10. [57] B. Gallien, M. Albaric, T. Duffar, K. Kakimoto, M. M’Hamdi, Study on the usage of a commercial software (Comsol-Multiphysics®) for dislocation multiplication model, J. Cryst. Growth ISSN (2016) ISSN 0022-0248. [58] B. Gao, K. Kakimoto, Dislocation-density-based modeling of the plastic behavior of 4H-SiC single crystals using the Alexander-Haasen model, J. Cryst. Growth 386 (2014) 215–219. [59] J. Steiner, M. Roder, B.D. Nguyen, S. Sandfeld, A. Danilewsky, P.J. Wellmann, Analysis of the basal plane dislocation density and thermomechanical stress during 100 mm pvt growth of 4h-sic, Materials 12 (13) (2019) ISSN 1996–1944..

10