- Email: [email protected]
Cellular-automata-based simulation of anisotropic crystal growth
Journal of Crystal Growth 220 (2000) 471–479
Cellular-automata-based simulation of anisotropic crystal growth M. Chahoud1, D. Fehly, H.-H. Wehmann*, A. Schlachetzki Institut fu¨r Halbleitertechnik, Technische Universita¨t Carolo-Wilhelmina zu Braunschweig, Postfach 3329, D-38023 Braunschweig, Germany Received 22 December 1998; accepted 5 October 2000 Communicated by J.J. Derby
Abstract Extending the simulation of anisotropic etching, a cellular-automata-based simulator is applied to anisotropic crystal growth. This simulator takes advantage of the equivalence between dissolution and growth of crystals. Metalorganic vapour-phase epitaxial growth experiments were performed on patterned (1 0 0)-oriented InP substrates with very deep V-shaped grooves with {1 1 1}A sidewalls to determine the relevant growth rates of InGaAs and InP. The capability of the simulation method is demonstrated by quantitative comparison of simulated and experimental results. In addition, the versatility of the model is shown with area-selective growth. # 2000 Elsevier Science B.V. All rights reserved. PACS: 07.05.T; 81.15.K; 81.05.D; 68.55 Keywords: Simulation; Cellular automata; Metalorganic vapour-phase epitaxy; Area-selective growth; Patterned substrate
1. Introduction Crystal growth on patterned or partly masked substrates plays an important role in the fabrication of semiconductor devices like waveguides and quantum-wire lasers [1–3]. The result of the growth process depends on many factors such as growth method, conditions and duration, crystal and surface structure, etc. Thus, the prediction of the resulting patterns is sometimes very difficult and a reliable growth simulator would minimize *Corresponding author. Tel: +49-531-391-3760, fax: +49531-391-5844. E-mail address: [email protected] (H.-H. Wehmann). 1 Present address: SMTP, AECS, P.O. Box 6091, Damascus, Syria.
the technological effort in fabricating new structures. In several studies the growth behaviour of InP-based semiconductors on planar, partially masked or structured substrates was investigated [3–10]. Two factors were found to be important for the growth rate: *
*
the orientation of the crystal face, i.e. the anisotropy of growth rates; the diffusion phenomena in the gas phase or at the surface of the crystal.
The growth parameters, mainly temperature and pressure determine which one dominates [9,10]. Even the applied growth method influences the surface shape. In near-equilibrium vapourphase epitaxy the interaction of growth species
0022-0248/00/$ - see front matter # 2000 Elsevier Science B.V. All rights reserved. PII: S 0 0 2 2 - 0 2 4 8 ( 0 0 ) 0 0 9 0 2 - 7
472
M. Chahoud et al. / Journal of Crystal Growth 220 (2000) 471–479
with the surface plays an important role. Thus, the growth rate is considerably influenced by dislocations and the orientation of surface steps. Whereas molecular beam epitaxy is more distant from equilibrium and the surface structure has less influence. Instead, surface migration determines the growth [11]. Beyond this, other parameters like lattice mismatch lead to growth instabilities with surface roughening and island growth which degrade the performance especially of quantumfilm devices. However, none of the currently used simulation programmes cover the whole field from ‘‘microscopic’’ to ‘‘macroscopic’’ features [12]. *
*
On the atomistic scale, in principle, there are two different categories – an analysis-based method and kinetic Monte Carlo simulations – which both may involve diffusion phenomena. The analytical approach uses homogeneous rate equations and continuum equations of motion. Since for practical applications the required number of parameters grows quickly, this method is primarily used for performing stability analyses [13] and identifying asymptotic regimes (scaling) [14]. On the other hand, the Monte Carlo approach is usually based on the length and time scales of single atoms and thus is only feasible for small details of complete device structures [15]. The level-set method is a combination of both, in a sense that the crystal is treated to be discrete in growth direction and continuous in lateral direction. It has the advantage of arbitrarily large lateral length scales [12]. However, since overhangs and undercuts are excluded, it is not applicable for areaselective growth of III/V-semiconductors [4]. Concerning ‘‘macroscopic’’ structures, geometrical models are preferred which simulate the motion of crystal facets during growth [16]. These may even include surface migration of adatoms [17]. With desirable three-dimensional structures, besides the increasing complexity in following all emerging planes and their intersection lines, a difficulty appears, how to handle saddle points where convex and concave edges meet. This problem has only been solved for crystal dissolution [18]. However, the implementation of diffusive transport in a three-
dimensional simulator is more difficult. Thus, only two-dimensional cases have been tackled yet [17,19]. In this paper, we concentrate on three-dimensional and ‘‘macroscopic’’ problems taking advantage of the equivalence of the dissolution of crystals during etching and crystal growth [20,21] in order to develop a new simulator for orientation-dependent growth which is relevant for patterned substrates and/or area-selective growth. This simulator is based on a deterministic cellularautomata etching model [22,23] which simulates the motion of facets during anisotropic etching as well as their appearance and disappearance in three dimensions of arbitrary systems. The only restriction is the neglection of surface migration of adatoms, since atomistic effects are globally taken into account by anisotropic growth rates. The neglect of diffusion can be done if one or more of the following assumptions are satisfied: (1) near equilibrium growth, (2) low temperatures, (3) fully faceted Wulff shape, (4) negligible lattice mismatch, (5) large structures as compared to migration lengths. The only experimental data needed for the simulation are the maxima and minima in the growth-rate polar diagram, which represents the growth-rate dependence on the crystallographic orientation. These minima and maxima dominate in convex and concave structures, respectively [20]. In order to examine the quality of this simulation method, a quantitative comparison between experimental and simulated structures is performed. The capability of this model to simulate the anisotropy in selective-area growth is also demonstrated.
2. Simulation model Because of the equivalence of etching concave and growing convex structures on the one hand and etching convex and growing concave structures on the other hand [20,23], we use the etching simulator developed in Refs. [22,23] to simulate the anisotropic growth after inverting the structures on which the growth process will be performed. The etching procedure is based on a
M. Chahoud et al. / Journal of Crystal Growth 220 (2000) 471–479
473
deterministic cellular-automata model. In this model [22,23]: 1. the crystal is represented by cells in a cubic primitive lattice; 2. the etching progress is simulated stepwise; 3. in each etching step cells within the surface are examined concerning their belonging to specific crystallographic planes; 4. a cell belongs to a crystallographic plane if all its neighbouring cells in this plane are present; 5. according to the corresponding etching rates it is decided whether the cell under test is immune against etching attack or not; 6. if it is not immune, the cell is removed from the solid; 7. higher etching rates of certain planes are accounted for by treating these planes more frequently. Although the simulator is applicable to threedimensional structures, in the following we illustrate this method by the two-dimensional example of Fig. 1. This figure shows convex ( , ) and concave ( , ) corners. For simplicity we consider only the {1 0} and {1 1} planes. The third Miller index is omitted since we deal with a twodimensional problem. In this example we assume the growth rate of {1 0} planes as d/3, while that of the {1 1} planes is d. At first, we symbolically invert the structure of Fig. 1a, i.e. white becomes dark and vice versa. We again stress that this inversion is just a geometrical procedure in order to employ the etching algorithm which has proven its success. In this particular example we substitute the growth of convex structures, with the slowest planes dominant, and of concave structures, with the fastest planes dominant, by etching of concave structures (slowest planes) and of convex structures (fastest planes), respectively. We obtain the pattern of Fig. 1b, now with concave corners instead of convex ones at and and with convex corners at and . Because of the equivalence mentioned above this structure has to be etched with etching rates for the {1 0} and {1 1} planes which are equal to their growth rates (d/3 for the {1 0} and d for the {1 1} planes). After etching of the structure in Fig. 1b the {1 0} planes have moved relatively little, whereas the
Fig. 1. Schematic illustration of the simulation algorithm: (a) original structure of crystal; (b) inverted structure for etching; (c) reinverted pattern, showing grown structure (bold lines).
{1 1} planes are moved by the distance d after one unit time. At the concave corners and the slowest planes, i. e. the {1 0} planes, are dominant [20,21]. Therefore, their shape remains unchanged. and the However, at the convex corners fastest planes, the {1 1} planes in this case, determine the final shape. Thus, the lightly shaded areas are dissolved. In Fig. 1c the etched structure is reinverted resulting in the final shape of the grown structure as marked by the bold lines. The extension of the simulation method to the threedimensional case is straightforward. Since the simulator is cellular-based the additional expenditure towards the third dimension is much less than in case of geometrical simulators which have to keep track of the motion of every present plane as well as of emerging facets and their intersection lines. In addition to that, the purely deterministic treatment in connection with the neglect of surface migration leads to a limited number of cells which have to be considered in contrast to Monte Carlo methods where a huge number of particles has to be treated in order to get meaningful statistics.
474
M. Chahoud et al. / Journal of Crystal Growth 220 (2000) 471–479
3. Experimental In order to examine the quality of the simulation of anisotropic growth processes on patterned substrates the growth rates of InGaAs and InP in h1 1 1i and h1 0 0i directions were determined experimentally. We used semi-insulating (1 0 0)oriented InP substrates with V-grooves etched through a titanium mask with a two-step etching process in order to obtain very smooth sidewalls [24]. Their width as well as their distance were varied between 5 and 200 mm. For the first etching step we used 47% HBr for 10 min followed by the second for 5 s in HBr:K2Cr2O7 (3:1) which smoothed out the residual roughness of the {1 1 1}A planes making them suitable for the metalorganic vapour-phase epitaxial growth (MOVPE) which was carried out at a temperature of 6008C and a total pressure of 100 hPa. The layer stack grown consisted of a 5 nm thick InP buffer layer, followed by 1 mm nominally lattice-matched InGaAs, 0.5 mm InP, 0.1 mm again lattice-matched InGaAs and finally 0.1 mm InP. We chose the growth conditions such that we expect a composition In0.53Ga0.47As, which we denote as InGaAs. The given thicknesses as well as the intended ternary lattice-matched composition are nominal values obtained on a flat (1 0 0)InP substrate. The actual orientation-dependent layer thicknesses were determined by scanning electron microscopy (SEM) on the cleavage plane perpendicular to the V-grooves as shown schematically in the central part of Fig. 2. The material contrast was established by selective etching of InGaAs for 1 min in citric acid:H2O2 (7:1) [25]. The SEM images of Fig. 2a–d were taken at the most interesting positions of the cleaved sample. The growth rates of the {1 0 0} and {1 1 1}A planes were determined at positions which are in the middle of the 200 mm wide mesas (Fig. 2b) between adjacent V-grooves and of the mesa ramps of the 100 mm wide grooves (Fig. 2d). At these positions the lateral diffusion from the neighbouring crystallographic planes is very weak, because these planes are very far away. In the tip of the V–grooves {3 1 1}A planes appear. Because of their small width (only about 1 mm), their growth rates are not primarily determined by the crystal anisotropy but
are also strongly influenced by diffusion from the neighbouring {1 1 1}A planes. The resulting ‘‘combined’’ growth rates are displayed in Table 1. An interesting point should be mentioned concerning the thick ternary layer (dark) in the SEM images of Fig. 2a and c. Clearly visible is a straight line between the lower edges to the underlying InP buffer and the upper edges to the second InP layer separating two areas of different brightness. For clarity we marked the beginning and the end of the respective lines by white crosses. The lines delineate the path of the moving edge during growth and separate the growing (1 0 0) and {1 1 1}A layers. Since the etch rate of the materialselective etchant used in the present case is sensitive to the composition of InGaAs [26] we attribute the contrast in brightness to a different Ga-content in the ternary layer grown on {1 1 1}A planes as opposed to the (1 0 0) surface. This is in accordance with findings on planar substrates of different orientation [27] as well as with a similar structure [28] as ours.
4. Results Starting with the experimentally determined growth rates of Table 1 we compare in Fig. 3 the experiment (part a) with our simulation (part b) of a cleaved mesa structure with the same layer sequence as in Fig. 2. The computing time of this two-dimensional pattern with cell dimensions of (20 nm)2 was about 3 min on a PC (Pentium 120 MHz). In order to distinguish between the different materials and compositions we use different shadings for InP cells (light gray) and InGaAs cells grown on (1 0 0) (dark gray) and {1 1 1}A planes (medium gray). For a quantitative comparison the simulated and experimental values of the thickness of the thicker binary and ternary layers (InP and InGaAs) in both orientations are listed in Table 2. Additionally, the angle a between the [0 1 1]-direction and the borderline between the InGaAs areas of different composition is given. With the sole exception of the (1 0 0) InGaAs layer thickness the agreement between simulation and experiment is good, including the borderline which could be simulated for the first time.
M. Chahoud et al. / Journal of Crystal Growth 220 (2000) 471–479
The most propable reason for the deviation in (1 0 0) layer thickness is the lateral diffusion from the {1 1 1}A planes. This diffusion process is caused by the lower {1 1 1}A growth rates and thus a higher concentration of mobile growth species on the sidewalls. The amount of deviation due to diffusion strongly depends on the geometry and is negligible for narrow V-grooves [29]. This effect does not appear with InP because its anisotropy in the growth rates under the applied
475
parameters is much less than that for InGaAs (cf. Table 1). As with our etching simulator, diffusion processes are not included in our purely deterministic model.
5. Extension of model The extension of the simulation model to area-selective growth is straightforward. Fig. 4
Fig. 2. SEM images at different positions of the grown and cleaved sample.
476
M. Chahoud et al. / Journal of Crystal Growth 220 (2000) 471–479
demonstrates this possibility with a fictitious set of growth rates. The growth of a diamond lattice structure (like Si) through a square window of (90 mm)2 is modeled in part a. The mask (white lines) is aligned along the [0 1 0] direction. The growth rates are assumed to correspond to the etching rates of silicon in KOH–isopropanol [30]. The simulated growth pattern reflects the four-fold symmetry of the diamond lattice. Taking the same mask as starting point, Fig. 4b shows the simulation for a zincblende lattice based on growth rates corresponding to etching rates for InP in HBr given in Ref. [22]. The result clearly exhibits the reduced symmetry as compared to the upper picture due to different growth rates in the [0 1 1]
and [0 1 1] directions. Although no experimental proof is available, Fig. 4 demonstrates first the ability of our model for three-dimensional simulations and second its versatility for arbitrary crystal types.
Table 1 Growth rates of InGaAs and InP as determined from Fig. 2 Growth rate InGaAs (mm/h) InP (mm/h)
{1 0 0} plane 2 0.1 1.06 0.04
{1 1 1}A plane
{3 1 1}A plane
0.66 0.1
1.6 0.2
0.94 0.08
1.12 0.14
Fig. 3. Experimental (a) and simulation results (b) of a grown mesa structure. (a) SEM microphotograph of the cleaved and etched (0 1 1) plane.
477
M. Chahoud et al. / Journal of Crystal Growth 220 (2000) 471–479
Table 2 Quantitative comparison between simulation and experiment of thicknesses of the two layers (indices b and t binary and ternary, respectively) grown first and the angle a in Fig. 3 ð1 0 0Þ
db Experiment Simulation
ðmmÞ
0.51 0.09 0.5
ð1 0 0Þ
dt
ðmmÞ
1.43 0.12 1
ð1 1 1Þ
db
ðmmÞ
0.49 0.1 0.45
ð1 1 1Þ
dt
ðmmÞ
0.39 0.08 0.33
a (8) 69 75
Fig. 4. Simulation of three-dimensional area-selective growth through a square window with fictitious growth rates: (a) diamond structure; (b) zincblende structure. The contours on the right margin give the sections AA’ and BB’, respectively.
For a quantitative comparison with experiments we restrict ourselves again to two-dimensional structures. Fig. 5a shows an SEM image of the (0 1 1) cleavage plane of an area-selectively grown, lattice-mismatched layer stack, whereas in part b the respective simulation is depicted. The maskless selective growth process was performed similar to [31] on a 1.5 mm GaAs layer overgrown with a lowtemperature 40 nm thick InP nucleation layer,
both on a (1 0 0) Si substrate. This double layer was partly etched down to the substrate by applying the same stripe mask as described previously. Thus, the unetched areas serve as nuclei for the following selective growth process which was performed by MOVPE at 6408C and 20 hPa. In this second epitaxial step the following layer sequence was grown with the nominal, i.e. on planar (1 0 0) substrates expected, layer
478
M. Chahoud et al. / Journal of Crystal Growth 220 (2000) 471–479
Fig. 5. Experimental (a) and simulation (b) results of a two-dimensional area-selective growth process.
Table 3 Quantitative comparison between simulation and experiment of Fig. 5
Experiment Simulation dexp =dsim :
ð1 0 0Þ
ðmmÞ
a (8)
b (8)
d1
103.1 4 100.2 }
167.2 2 177.1 }
4.1 0.2 1.88 2.18
thicknesses: InP (1.9 mm), InP–InGaAs multiquantum well (20 periods of 7 nm InP/7 nm InGaAs), InP (1.2 mm), InGaAsP (0.2 mm), and finally InP (0.9 mm). For visualization of the layer sequence we adopted material-selective etching by H2SO4 :H2O2 : H2O (5 : 1 : 1) for 20 s. For the simulation the {1 0 0} growth rates of InP and InGaAs were taken from Ref. [4]. The respective value of InGaAsP was determined on a simultaneously grown layer on a planar InP substrate and equals 1.9 mm/h. For the {1 1 1}A and B planes there is no information about the growth rate available at the applied process parameters. Assuming the ratio between the rates on {1 1 1}A and {1 0 0} planes to be constant even if temperature is varied [32], we calculated the {1 1 1}A growth rate using the ratio taken from
ð1 1 1ÞA
d1
ðmmÞ
3.2 0.15 1.4 2.3
ð1 1 1ÞB
d1
ðmmÞ
2.88 0.2 1.58 1.8
Table 1 and the {1 0 0} rate given in Ref. [4]. A similar procedure leads to the growth rate of {1 1 1}B planes. However, in this case the growth ratio between both {1 1 1} planes was taken from Ref. [4] where smaller structures where examined and thus it is more strongly influenced by diffusion effects. Comparing experiment and simulation in Fig. 5, in both pictures the {1 1 1}A and B planes are clearly visible and the qualitative agreement is good. For a quantitative comparison we collected in Table 3 the layer thicknesses taken on the three planes marked in part b as well as two angles as given in part a. The agreement for a and b is excellent, whereas the thicknesses between simulation and experiment differ by about a factor of 2. This deviation again is caused by diffusion
M. Chahoud et al. / Journal of Crystal Growth 220 (2000) 471–479
processes from the neighbouring Si surface on which no nucleation takes place at the given growth conditions. 6. Conclusion A model for three-dimensional anisotropic growth simulations is developed on the basis of a cellular-automata-based etching simulator which can be employed when growing on patterned substrates or on selected areas. The deterministic approach allows the efficient simulation also in cases the growth-rate anisotropy is only small. A quantitative comparison between experimental and simulation results shows that the anisotropy of growth is successfully simulated. Even the border line between different compositions of InGaAs originating from growth on different crystal planes could be simulated for the first time. However, deviations between simulation and experiment occur, when diffusion phenomena are involved. Acknowledgements The authors wish to thank D. Wu¨llner, P. Bo¨nsch, M. Karsten, and D. Ru¨mmler for helpful discussions and technical assitance. The valuable hints of the referees are acknowledged. M. C. thanks the Energy Research Group of the University of Damaskus for financial support. The assistance by the Deutsche Forschungsgemeinschaft and the Volkswagen Stiftung for parts of the experiments is gratefully acknowledged.
References [1] S. Tsukamoto, Y. Nagamune, M. Nishioka, Y. Arakawa, J. Appl. Phys. 71 (1992) 533. [2] M.L. Osowski, R. Panepucci, D.A. Turnbull, S.Q. Gu, A.M. Jones, S.G. Bishop, I. Adesida, J.J. Coleman, Appl. Phys. Lett. 68 (1996) 1087. [3] M. Kappelt, M. Grundmann, A. Krost, V. Tu¨rck, D. Bimberg, Appl. Phys. Lett. 68 (1996) 3596. [4] G. Zwinge, H.-H. Wehmann, A. Schlachetzki, C.C. Hsu, J. Appl. Phys. 74 (1993) 5516.
479
[5] T. Fujii, M. Ekawa, S. Yamazaki, J. Crystal Growth 156 (1995) 59. [6] T. Fujii, M. Ekawa, S. Yamazaki, J. Crystal Growth 146 (1995) 475. [7] M.J. Anders, M.M.G. Bongers, P.L. Bastos, L.J. Giling, J. Crystal Growth 154 (1995) 240. [8] K. Nakai, T. Sanada, S. Yamakoshi, J. Crystal Growth 93 (1988) 248. [9] T. Fujii, M. Ekawa, J. Appl. Phys. 78 (1995) 5373. [10] D.H. Reep, S.K. Ghandhi, J. Electrochem. Soc. 130 (1983) 675. [11] J. van Wingerden, R.H. van Aken, Y.A. Wiechers, P.M.L.O. Scholte, F. Tuinstra, Phys. Rev. B 57 (1998) 7252. [12] M.F. Gyure, C. Ratsch, B. Merriman, R.E. Caflisch, S. Osher, J.J. Zinck, D.D. Vvedensky, Phys. Rev. E 58 (1998) R6927. [13] L. Golubovic´, Phys. Rev. Lett. 78 (1997) 90. [14] T. Sun, M. Plischke, Phys. Rev. E 50 (1994) 3370. [15] H. Huang, G.H. Gilmer, T. Dı´ az de la Rubia, J. Appl. Phys. 84 (1998) 3636. [16] A. Bakin, G. Zwinge, R. Klockenbrink, H.-H. Wehmann, A. Schlachetzki, J. Appl. Phys. 76 (1994) 4906. [17] M. Ohtsuka, A. Suzuki, J. Appl. Phys. 73 (1993) 7358. [18] C.H. Se´quin, Sensors and Actuators A 34 (1992) 225. [19] W.C. Carter, A.R. Roosen, J.W. Cahn, J.E. Taylor, Acta Metall. Mater. 43 (1995) 4309. [20] R. Lacmann, W. Franke, R. Heimann, J. Crystal Growth 26 (1974) 107. [21] B.W. Batterman, J. Appl. Phys. 28 (1957) 1236. [22] M. Chahoud, H.-H. Wehmann, A. Schlachetzki, Sensors and Actuators A 63 (1997) 141. [23] M. Chahoud, H.-H. Wehmann, A. Schlachetzki, Sensors and Actuators A 69 (1998) 251. [24] P. Bo¨nsch, D. Wu¨llner, T. Schrimpf, A. Schlachetzki, R. Lacmann, J. Electrochem. Soc. 145 (1998) 1273. [25] D. Wu¨llner, A. Schlachetzki, P. Bo¨nsch, H.-H. Wehmann, T. Schrimpf, R. Lacmann, S. Kipp, Mater. Sci. Eng. B51 (1998) 178. [26] G.C. DeSalvo, W.F. Tseng, J. Comas, J. Electrochem. Soc. 139 (1992) 831. [27] O. Kayser, J. Crystal Growth 107 (1991) 989. [28] M. Kappelt, V. Tu¨rck, O. Stier, D. Bimberg, H. Cerva, D. Stenkamp, P. Veit, T. Hempel, J. Christen, Proceedings of the European Workshop MOVPE VII, Berlin (June 811, 1997) paper E5. [29] D. Wu¨llner, M. Chahoud, T. Schrimpf, P. Bo¨nsch, H.-H. Wehmann, A. Schlachetzki, J. Appl. Phys. 85 (1999) 249. [30] I. Barycka, I. Zubel, Sensors and Actuators A 48 (1995) 229. [31] G.-P. Tang, E. Peiner, H.-H. Wehmann, A. Lubnow, G. Zwinge, A. Schlachetzki, J. Hergeth, J. Appl. Phys. 72 (1992) 4366. [32] G.K. Mayer, H.L. Offereins, H. Sandmaier, K. Ku¨hl, J. Electrochem. Soc. 137 (1990) 3947.
Cellular-automata-based simulation of anisotropic crystal growth M. Chahoud1, D. Fehly, H.-H. Wehmann*, A. Schlachetzki Institut fu¨r Halbleitertechnik, Technische Universita¨t Carolo-Wilhelmina zu Braunschweig, Postfach 3329, D-38023 Braunschweig, Germany Received 22 December 1998; accepted 5 October 2000 Communicated by J.J. Derby
Abstract Extending the simulation of anisotropic etching, a cellular-automata-based simulator is applied to anisotropic crystal growth. This simulator takes advantage of the equivalence between dissolution and growth of crystals. Metalorganic vapour-phase epitaxial growth experiments were performed on patterned (1 0 0)-oriented InP substrates with very deep V-shaped grooves with {1 1 1}A sidewalls to determine the relevant growth rates of InGaAs and InP. The capability of the simulation method is demonstrated by quantitative comparison of simulated and experimental results. In addition, the versatility of the model is shown with area-selective growth. # 2000 Elsevier Science B.V. All rights reserved. PACS: 07.05.T; 81.15.K; 81.05.D; 68.55 Keywords: Simulation; Cellular automata; Metalorganic vapour-phase epitaxy; Area-selective growth; Patterned substrate
1. Introduction Crystal growth on patterned or partly masked substrates plays an important role in the fabrication of semiconductor devices like waveguides and quantum-wire lasers [1–3]. The result of the growth process depends on many factors such as growth method, conditions and duration, crystal and surface structure, etc. Thus, the prediction of the resulting patterns is sometimes very difficult and a reliable growth simulator would minimize *Corresponding author. Tel: +49-531-391-3760, fax: +49531-391-5844. E-mail address: [email protected] (H.-H. Wehmann). 1 Present address: SMTP, AECS, P.O. Box 6091, Damascus, Syria.
the technological effort in fabricating new structures. In several studies the growth behaviour of InP-based semiconductors on planar, partially masked or structured substrates was investigated [3–10]. Two factors were found to be important for the growth rate: *
*
the orientation of the crystal face, i.e. the anisotropy of growth rates; the diffusion phenomena in the gas phase or at the surface of the crystal.
The growth parameters, mainly temperature and pressure determine which one dominates [9,10]. Even the applied growth method influences the surface shape. In near-equilibrium vapourphase epitaxy the interaction of growth species
0022-0248/00/$ - see front matter # 2000 Elsevier Science B.V. All rights reserved. PII: S 0 0 2 2 - 0 2 4 8 ( 0 0 ) 0 0 9 0 2 - 7
472
M. Chahoud et al. / Journal of Crystal Growth 220 (2000) 471–479
with the surface plays an important role. Thus, the growth rate is considerably influenced by dislocations and the orientation of surface steps. Whereas molecular beam epitaxy is more distant from equilibrium and the surface structure has less influence. Instead, surface migration determines the growth [11]. Beyond this, other parameters like lattice mismatch lead to growth instabilities with surface roughening and island growth which degrade the performance especially of quantumfilm devices. However, none of the currently used simulation programmes cover the whole field from ‘‘microscopic’’ to ‘‘macroscopic’’ features [12]. *
*
On the atomistic scale, in principle, there are two different categories – an analysis-based method and kinetic Monte Carlo simulations – which both may involve diffusion phenomena. The analytical approach uses homogeneous rate equations and continuum equations of motion. Since for practical applications the required number of parameters grows quickly, this method is primarily used for performing stability analyses [13] and identifying asymptotic regimes (scaling) [14]. On the other hand, the Monte Carlo approach is usually based on the length and time scales of single atoms and thus is only feasible for small details of complete device structures [15]. The level-set method is a combination of both, in a sense that the crystal is treated to be discrete in growth direction and continuous in lateral direction. It has the advantage of arbitrarily large lateral length scales [12]. However, since overhangs and undercuts are excluded, it is not applicable for areaselective growth of III/V-semiconductors [4]. Concerning ‘‘macroscopic’’ structures, geometrical models are preferred which simulate the motion of crystal facets during growth [16]. These may even include surface migration of adatoms [17]. With desirable three-dimensional structures, besides the increasing complexity in following all emerging planes and their intersection lines, a difficulty appears, how to handle saddle points where convex and concave edges meet. This problem has only been solved for crystal dissolution [18]. However, the implementation of diffusive transport in a three-
dimensional simulator is more difficult. Thus, only two-dimensional cases have been tackled yet [17,19]. In this paper, we concentrate on three-dimensional and ‘‘macroscopic’’ problems taking advantage of the equivalence of the dissolution of crystals during etching and crystal growth [20,21] in order to develop a new simulator for orientation-dependent growth which is relevant for patterned substrates and/or area-selective growth. This simulator is based on a deterministic cellularautomata etching model [22,23] which simulates the motion of facets during anisotropic etching as well as their appearance and disappearance in three dimensions of arbitrary systems. The only restriction is the neglection of surface migration of adatoms, since atomistic effects are globally taken into account by anisotropic growth rates. The neglect of diffusion can be done if one or more of the following assumptions are satisfied: (1) near equilibrium growth, (2) low temperatures, (3) fully faceted Wulff shape, (4) negligible lattice mismatch, (5) large structures as compared to migration lengths. The only experimental data needed for the simulation are the maxima and minima in the growth-rate polar diagram, which represents the growth-rate dependence on the crystallographic orientation. These minima and maxima dominate in convex and concave structures, respectively [20]. In order to examine the quality of this simulation method, a quantitative comparison between experimental and simulated structures is performed. The capability of this model to simulate the anisotropy in selective-area growth is also demonstrated.
2. Simulation model Because of the equivalence of etching concave and growing convex structures on the one hand and etching convex and growing concave structures on the other hand [20,23], we use the etching simulator developed in Refs. [22,23] to simulate the anisotropic growth after inverting the structures on which the growth process will be performed. The etching procedure is based on a
M. Chahoud et al. / Journal of Crystal Growth 220 (2000) 471–479
473
deterministic cellular-automata model. In this model [22,23]: 1. the crystal is represented by cells in a cubic primitive lattice; 2. the etching progress is simulated stepwise; 3. in each etching step cells within the surface are examined concerning their belonging to specific crystallographic planes; 4. a cell belongs to a crystallographic plane if all its neighbouring cells in this plane are present; 5. according to the corresponding etching rates it is decided whether the cell under test is immune against etching attack or not; 6. if it is not immune, the cell is removed from the solid; 7. higher etching rates of certain planes are accounted for by treating these planes more frequently. Although the simulator is applicable to threedimensional structures, in the following we illustrate this method by the two-dimensional example of Fig. 1. This figure shows convex ( , ) and concave ( , ) corners. For simplicity we consider only the {1 0} and {1 1} planes. The third Miller index is omitted since we deal with a twodimensional problem. In this example we assume the growth rate of {1 0} planes as d/3, while that of the {1 1} planes is d. At first, we symbolically invert the structure of Fig. 1a, i.e. white becomes dark and vice versa. We again stress that this inversion is just a geometrical procedure in order to employ the etching algorithm which has proven its success. In this particular example we substitute the growth of convex structures, with the slowest planes dominant, and of concave structures, with the fastest planes dominant, by etching of concave structures (slowest planes) and of convex structures (fastest planes), respectively. We obtain the pattern of Fig. 1b, now with concave corners instead of convex ones at and and with convex corners at and . Because of the equivalence mentioned above this structure has to be etched with etching rates for the {1 0} and {1 1} planes which are equal to their growth rates (d/3 for the {1 0} and d for the {1 1} planes). After etching of the structure in Fig. 1b the {1 0} planes have moved relatively little, whereas the
Fig. 1. Schematic illustration of the simulation algorithm: (a) original structure of crystal; (b) inverted structure for etching; (c) reinverted pattern, showing grown structure (bold lines).
{1 1} planes are moved by the distance d after one unit time. At the concave corners and the slowest planes, i. e. the {1 0} planes, are dominant [20,21]. Therefore, their shape remains unchanged. and the However, at the convex corners fastest planes, the {1 1} planes in this case, determine the final shape. Thus, the lightly shaded areas are dissolved. In Fig. 1c the etched structure is reinverted resulting in the final shape of the grown structure as marked by the bold lines. The extension of the simulation method to the threedimensional case is straightforward. Since the simulator is cellular-based the additional expenditure towards the third dimension is much less than in case of geometrical simulators which have to keep track of the motion of every present plane as well as of emerging facets and their intersection lines. In addition to that, the purely deterministic treatment in connection with the neglect of surface migration leads to a limited number of cells which have to be considered in contrast to Monte Carlo methods where a huge number of particles has to be treated in order to get meaningful statistics.
474
M. Chahoud et al. / Journal of Crystal Growth 220 (2000) 471–479
3. Experimental In order to examine the quality of the simulation of anisotropic growth processes on patterned substrates the growth rates of InGaAs and InP in h1 1 1i and h1 0 0i directions were determined experimentally. We used semi-insulating (1 0 0)oriented InP substrates with V-grooves etched through a titanium mask with a two-step etching process in order to obtain very smooth sidewalls [24]. Their width as well as their distance were varied between 5 and 200 mm. For the first etching step we used 47% HBr for 10 min followed by the second for 5 s in HBr:K2Cr2O7 (3:1) which smoothed out the residual roughness of the {1 1 1}A planes making them suitable for the metalorganic vapour-phase epitaxial growth (MOVPE) which was carried out at a temperature of 6008C and a total pressure of 100 hPa. The layer stack grown consisted of a 5 nm thick InP buffer layer, followed by 1 mm nominally lattice-matched InGaAs, 0.5 mm InP, 0.1 mm again lattice-matched InGaAs and finally 0.1 mm InP. We chose the growth conditions such that we expect a composition In0.53Ga0.47As, which we denote as InGaAs. The given thicknesses as well as the intended ternary lattice-matched composition are nominal values obtained on a flat (1 0 0)InP substrate. The actual orientation-dependent layer thicknesses were determined by scanning electron microscopy (SEM) on the cleavage plane perpendicular to the V-grooves as shown schematically in the central part of Fig. 2. The material contrast was established by selective etching of InGaAs for 1 min in citric acid:H2O2 (7:1) [25]. The SEM images of Fig. 2a–d were taken at the most interesting positions of the cleaved sample. The growth rates of the {1 0 0} and {1 1 1}A planes were determined at positions which are in the middle of the 200 mm wide mesas (Fig. 2b) between adjacent V-grooves and of the mesa ramps of the 100 mm wide grooves (Fig. 2d). At these positions the lateral diffusion from the neighbouring crystallographic planes is very weak, because these planes are very far away. In the tip of the V–grooves {3 1 1}A planes appear. Because of their small width (only about 1 mm), their growth rates are not primarily determined by the crystal anisotropy but
are also strongly influenced by diffusion from the neighbouring {1 1 1}A planes. The resulting ‘‘combined’’ growth rates are displayed in Table 1. An interesting point should be mentioned concerning the thick ternary layer (dark) in the SEM images of Fig. 2a and c. Clearly visible is a straight line between the lower edges to the underlying InP buffer and the upper edges to the second InP layer separating two areas of different brightness. For clarity we marked the beginning and the end of the respective lines by white crosses. The lines delineate the path of the moving edge during growth and separate the growing (1 0 0) and {1 1 1}A layers. Since the etch rate of the materialselective etchant used in the present case is sensitive to the composition of InGaAs [26] we attribute the contrast in brightness to a different Ga-content in the ternary layer grown on {1 1 1}A planes as opposed to the (1 0 0) surface. This is in accordance with findings on planar substrates of different orientation [27] as well as with a similar structure [28] as ours.
4. Results Starting with the experimentally determined growth rates of Table 1 we compare in Fig. 3 the experiment (part a) with our simulation (part b) of a cleaved mesa structure with the same layer sequence as in Fig. 2. The computing time of this two-dimensional pattern with cell dimensions of (20 nm)2 was about 3 min on a PC (Pentium 120 MHz). In order to distinguish between the different materials and compositions we use different shadings for InP cells (light gray) and InGaAs cells grown on (1 0 0) (dark gray) and {1 1 1}A planes (medium gray). For a quantitative comparison the simulated and experimental values of the thickness of the thicker binary and ternary layers (InP and InGaAs) in both orientations are listed in Table 2. Additionally, the angle a between the [0 1 1]-direction and the borderline between the InGaAs areas of different composition is given. With the sole exception of the (1 0 0) InGaAs layer thickness the agreement between simulation and experiment is good, including the borderline which could be simulated for the first time.
M. Chahoud et al. / Journal of Crystal Growth 220 (2000) 471–479
The most propable reason for the deviation in (1 0 0) layer thickness is the lateral diffusion from the {1 1 1}A planes. This diffusion process is caused by the lower {1 1 1}A growth rates and thus a higher concentration of mobile growth species on the sidewalls. The amount of deviation due to diffusion strongly depends on the geometry and is negligible for narrow V-grooves [29]. This effect does not appear with InP because its anisotropy in the growth rates under the applied
475
parameters is much less than that for InGaAs (cf. Table 1). As with our etching simulator, diffusion processes are not included in our purely deterministic model.
5. Extension of model The extension of the simulation model to area-selective growth is straightforward. Fig. 4
Fig. 2. SEM images at different positions of the grown and cleaved sample.
476
M. Chahoud et al. / Journal of Crystal Growth 220 (2000) 471–479
demonstrates this possibility with a fictitious set of growth rates. The growth of a diamond lattice structure (like Si) through a square window of (90 mm)2 is modeled in part a. The mask (white lines) is aligned along the [0 1 0] direction. The growth rates are assumed to correspond to the etching rates of silicon in KOH–isopropanol [30]. The simulated growth pattern reflects the four-fold symmetry of the diamond lattice. Taking the same mask as starting point, Fig. 4b shows the simulation for a zincblende lattice based on growth rates corresponding to etching rates for InP in HBr given in Ref. [22]. The result clearly exhibits the reduced symmetry as compared to the upper picture due to different growth rates in the [0 1 1]
and [0 1 1] directions. Although no experimental proof is available, Fig. 4 demonstrates first the ability of our model for three-dimensional simulations and second its versatility for arbitrary crystal types.
Table 1 Growth rates of InGaAs and InP as determined from Fig. 2 Growth rate InGaAs (mm/h) InP (mm/h)
{1 0 0} plane 2 0.1 1.06 0.04
{1 1 1}A plane
{3 1 1}A plane
0.66 0.1
1.6 0.2
0.94 0.08
1.12 0.14
Fig. 3. Experimental (a) and simulation results (b) of a grown mesa structure. (a) SEM microphotograph of the cleaved and etched (0 1 1) plane.
477
M. Chahoud et al. / Journal of Crystal Growth 220 (2000) 471–479
Table 2 Quantitative comparison between simulation and experiment of thicknesses of the two layers (indices b and t binary and ternary, respectively) grown first and the angle a in Fig. 3 ð1 0 0Þ
db Experiment Simulation
ðmmÞ
0.51 0.09 0.5
ð1 0 0Þ
dt
ðmmÞ
1.43 0.12 1
ð1 1 1Þ
db
ðmmÞ
0.49 0.1 0.45
ð1 1 1Þ
dt
ðmmÞ
0.39 0.08 0.33
a (8) 69 75
Fig. 4. Simulation of three-dimensional area-selective growth through a square window with fictitious growth rates: (a) diamond structure; (b) zincblende structure. The contours on the right margin give the sections AA’ and BB’, respectively.
For a quantitative comparison with experiments we restrict ourselves again to two-dimensional structures. Fig. 5a shows an SEM image of the (0 1 1) cleavage plane of an area-selectively grown, lattice-mismatched layer stack, whereas in part b the respective simulation is depicted. The maskless selective growth process was performed similar to [31] on a 1.5 mm GaAs layer overgrown with a lowtemperature 40 nm thick InP nucleation layer,
both on a (1 0 0) Si substrate. This double layer was partly etched down to the substrate by applying the same stripe mask as described previously. Thus, the unetched areas serve as nuclei for the following selective growth process which was performed by MOVPE at 6408C and 20 hPa. In this second epitaxial step the following layer sequence was grown with the nominal, i.e. on planar (1 0 0) substrates expected, layer
478
M. Chahoud et al. / Journal of Crystal Growth 220 (2000) 471–479
Fig. 5. Experimental (a) and simulation (b) results of a two-dimensional area-selective growth process.
Table 3 Quantitative comparison between simulation and experiment of Fig. 5
Experiment Simulation dexp =dsim :
ð1 0 0Þ
ðmmÞ
a (8)
b (8)
d1
103.1 4 100.2 }
167.2 2 177.1 }
4.1 0.2 1.88 2.18
thicknesses: InP (1.9 mm), InP–InGaAs multiquantum well (20 periods of 7 nm InP/7 nm InGaAs), InP (1.2 mm), InGaAsP (0.2 mm), and finally InP (0.9 mm). For visualization of the layer sequence we adopted material-selective etching by H2SO4 :H2O2 : H2O (5 : 1 : 1) for 20 s. For the simulation the {1 0 0} growth rates of InP and InGaAs were taken from Ref. [4]. The respective value of InGaAsP was determined on a simultaneously grown layer on a planar InP substrate and equals 1.9 mm/h. For the {1 1 1}A and B planes there is no information about the growth rate available at the applied process parameters. Assuming the ratio between the rates on {1 1 1}A and {1 0 0} planes to be constant even if temperature is varied [32], we calculated the {1 1 1}A growth rate using the ratio taken from
ð1 1 1ÞA
d1
ðmmÞ
3.2 0.15 1.4 2.3
ð1 1 1ÞB
d1
ðmmÞ
2.88 0.2 1.58 1.8
Table 1 and the {1 0 0} rate given in Ref. [4]. A similar procedure leads to the growth rate of {1 1 1}B planes. However, in this case the growth ratio between both {1 1 1} planes was taken from Ref. [4] where smaller structures where examined and thus it is more strongly influenced by diffusion effects. Comparing experiment and simulation in Fig. 5, in both pictures the {1 1 1}A and B planes are clearly visible and the qualitative agreement is good. For a quantitative comparison we collected in Table 3 the layer thicknesses taken on the three planes marked in part b as well as two angles as given in part a. The agreement for a and b is excellent, whereas the thicknesses between simulation and experiment differ by about a factor of 2. This deviation again is caused by diffusion
M. Chahoud et al. / Journal of Crystal Growth 220 (2000) 471–479
processes from the neighbouring Si surface on which no nucleation takes place at the given growth conditions. 6. Conclusion A model for three-dimensional anisotropic growth simulations is developed on the basis of a cellular-automata-based etching simulator which can be employed when growing on patterned substrates or on selected areas. The deterministic approach allows the efficient simulation also in cases the growth-rate anisotropy is only small. A quantitative comparison between experimental and simulation results shows that the anisotropy of growth is successfully simulated. Even the border line between different compositions of InGaAs originating from growth on different crystal planes could be simulated for the first time. However, deviations between simulation and experiment occur, when diffusion phenomena are involved. Acknowledgements The authors wish to thank D. Wu¨llner, P. Bo¨nsch, M. Karsten, and D. Ru¨mmler for helpful discussions and technical assitance. The valuable hints of the referees are acknowledged. M. C. thanks the Energy Research Group of the University of Damaskus for financial support. The assistance by the Deutsche Forschungsgemeinschaft and the Volkswagen Stiftung for parts of the experiments is gratefully acknowledged.
References [1] S. Tsukamoto, Y. Nagamune, M. Nishioka, Y. Arakawa, J. Appl. Phys. 71 (1992) 533. [2] M.L. Osowski, R. Panepucci, D.A. Turnbull, S.Q. Gu, A.M. Jones, S.G. Bishop, I. Adesida, J.J. Coleman, Appl. Phys. Lett. 68 (1996) 1087. [3] M. Kappelt, M. Grundmann, A. Krost, V. Tu¨rck, D. Bimberg, Appl. Phys. Lett. 68 (1996) 3596. [4] G. Zwinge, H.-H. Wehmann, A. Schlachetzki, C.C. Hsu, J. Appl. Phys. 74 (1993) 5516.
479
[5] T. Fujii, M. Ekawa, S. Yamazaki, J. Crystal Growth 156 (1995) 59. [6] T. Fujii, M. Ekawa, S. Yamazaki, J. Crystal Growth 146 (1995) 475. [7] M.J. Anders, M.M.G. Bongers, P.L. Bastos, L.J. Giling, J. Crystal Growth 154 (1995) 240. [8] K. Nakai, T. Sanada, S. Yamakoshi, J. Crystal Growth 93 (1988) 248. [9] T. Fujii, M. Ekawa, J. Appl. Phys. 78 (1995) 5373. [10] D.H. Reep, S.K. Ghandhi, J. Electrochem. Soc. 130 (1983) 675. [11] J. van Wingerden, R.H. van Aken, Y.A. Wiechers, P.M.L.O. Scholte, F. Tuinstra, Phys. Rev. B 57 (1998) 7252. [12] M.F. Gyure, C. Ratsch, B. Merriman, R.E. Caflisch, S. Osher, J.J. Zinck, D.D. Vvedensky, Phys. Rev. E 58 (1998) R6927. [13] L. Golubovic´, Phys. Rev. Lett. 78 (1997) 90. [14] T. Sun, M. Plischke, Phys. Rev. E 50 (1994) 3370. [15] H. Huang, G.H. Gilmer, T. Dı´ az de la Rubia, J. Appl. Phys. 84 (1998) 3636. [16] A. Bakin, G. Zwinge, R. Klockenbrink, H.-H. Wehmann, A. Schlachetzki, J. Appl. Phys. 76 (1994) 4906. [17] M. Ohtsuka, A. Suzuki, J. Appl. Phys. 73 (1993) 7358. [18] C.H. Se´quin, Sensors and Actuators A 34 (1992) 225. [19] W.C. Carter, A.R. Roosen, J.W. Cahn, J.E. Taylor, Acta Metall. Mater. 43 (1995) 4309. [20] R. Lacmann, W. Franke, R. Heimann, J. Crystal Growth 26 (1974) 107. [21] B.W. Batterman, J. Appl. Phys. 28 (1957) 1236. [22] M. Chahoud, H.-H. Wehmann, A. Schlachetzki, Sensors and Actuators A 63 (1997) 141. [23] M. Chahoud, H.-H. Wehmann, A. Schlachetzki, Sensors and Actuators A 69 (1998) 251. [24] P. Bo¨nsch, D. Wu¨llner, T. Schrimpf, A. Schlachetzki, R. Lacmann, J. Electrochem. Soc. 145 (1998) 1273. [25] D. Wu¨llner, A. Schlachetzki, P. Bo¨nsch, H.-H. Wehmann, T. Schrimpf, R. Lacmann, S. Kipp, Mater. Sci. Eng. B51 (1998) 178. [26] G.C. DeSalvo, W.F. Tseng, J. Comas, J. Electrochem. Soc. 139 (1992) 831. [27] O. Kayser, J. Crystal Growth 107 (1991) 989. [28] M. Kappelt, V. Tu¨rck, O. Stier, D. Bimberg, H. Cerva, D. Stenkamp, P. Veit, T. Hempel, J. Christen, Proceedings of the European Workshop MOVPE VII, Berlin (June 811, 1997) paper E5. [29] D. Wu¨llner, M. Chahoud, T. Schrimpf, P. Bo¨nsch, H.-H. Wehmann, A. Schlachetzki, J. Appl. Phys. 85 (1999) 249. [30] I. Barycka, I. Zubel, Sensors and Actuators A 48 (1995) 229. [31] G.-P. Tang, E. Peiner, H.-H. Wehmann, A. Lubnow, G. Zwinge, A. Schlachetzki, J. Hergeth, J. Appl. Phys. 72 (1992) 4366. [32] G.K. Mayer, H.L. Offereins, H. Sandmaier, K. Ku¨hl, J. Electrochem. Soc. 137 (1990) 3947.