Analytic many-body potential for InAs / GaAs surfaces and nanostructures : Formation energy of InAs quantum dots

A parametrization of the Abell–Tersoff potential for In, Ga, As, InAs, and GaAs is presented by using both experimental data and results from density-functional calculations as input. This parametrization is optimized for the description of structural and elastic properties of bulk In, Ga, As, InAs, and GaAs, as well as for the structure and energy of several reconstructed low-index GaAs and InAs surfaces. We demonstrate the transferability to GaAs and InAs high-index surfaces and compare the results to those obtained with previously published parametrizations. Furthermore, we demonstrate the applicability to epitaxial InAs/GaAs films by comparing the Poisson ratio and elastic energy for biaxial strain, as obtained numerically with our potential and analytically from continuum-elasticity theory. Limitations for the description of point defects and surface diffusion are pointed out. This parametrization enables us to perform atomically detailed studies of InAs/GaAs heterostructures. The formation energy of InAs quantum dots on GaAs 001 obtained from our atomistic approach is in good agreement with previous results from a hybrid approach.


I. INTRODUCTION
Contemporary semiconductor technology utilizes structures in which the mobility of the electronic carriers is restricted in one, two, or all three dimensions of space ͑quantum wells, wires, and dots, respectively͒.The typical length scale of the confining dimension͑s͒ can be as small as only a few atomic layers.For computational material science, this poses the challenge of atomistic modeling since the position of each atom may matter for the properties of these structures.At the same time, the strain fields associated with lattice-mismatched heterostructures can be rather long ranged such that many atoms need to be included in the simulations.Hence, a simulation technique that is accurate yet computationally not too expensive is sought for.In particular, such a technique should yield accurate atomic positions and strain fields since these are required as input for electronic structure calculations ͑using, e.g., effective mass, k • p theory, or tight-binding approaches͒ to calculate the electronic properties of the heterostructures.If the dimensions of the heterostructures are sufficiently large, continuumelasticity theory ͑CET͒ can be used to calculate strain fields, which are then used as input to continuum approaches to electronic structure, such as k • p theory.If an atomistic treatment is required, various types of valence force fields have been in use to determine relaxed atomic positions in strained heterostructures ͑see, e.g., Refs.1-3͒.It has been shown that contact with macroscopic elasticity theory can be made by a suitably defined strain tensor derived from the atomic coordinates. 4or self-assembled quantum dots or wires, computational modeling could help to elucidate the stability of different possible structures and to simulate elementary atomic processes during their growth.It is challenging to devise a suitable computationally inexpensive simulation technique for addressing these questions: free surfaces and other situations of low atomic coordination obviously play an important role in the fabrication of the nanostructures.However, commonly used valence force fields are parametrized only for bulk properties and cannot reliably describe the bonding at surfaces.Therefore, improved force fields are needed that include reliable microscopic information about surfaces.
In this work, we develop and test an analytical many-body potential that is particularly suited for the simulation of InAs and GaAs nanostructures.We are interested in the early stages of three-dimensional InAs islands grown on GaAs in the Stranski-Krastanov growth mode.After overgrowth with a suitable capping layer, these islands can be used as quantum dots.Several aspects of this system require atomistic modeling: for example, the critical thickness of the wetting layer for InAs/GaAs is less than two monolayers, which calls for an atomistic description of the wetting layer.Among the different types of potentials in use, we decided for the manybody potential of the Abell-Tersoff type, 5,6 as it is well suited for semiconductors with covalent bonding.We demand a highly accurate description of the lattice constants and elastic moduli of both GaAs and InAs, as well as a good description of the surface reconstructions and surface energies that occur under typical growth conditions for InAs/GaAs quantum dots.
Abell-Tersoff-type potentials depend on bond lengths and bond angles and hence access information about the atomic structure of a given material via the arrangements of pairs and triples of atoms.The dependence on these input data is usually parametrized for the ground state atomic structure of a chemical element or binary compound.It is not a priori clear that a particular parametrization gives reliable results when applied to new and different structures.For instance, in ternary compounds, additional triples of atoms not present in binary compounds occur.Moreover, significant deviations of the bond lengths and bond angles from those in bulk structures appear near point defects or at surfaces.It is not clear if a particular parametrization can reliably describe such situations that significantly deviate from the structures used for its derivation.Therefore, parametrizations employed in the literature to treat a specific simulation task should not be trans-ferred to a new modeling task without performing careful tests.Moreover, there are inherent limits to the performance of these potentials: since the functional form of Tersoff-Abell potentials offers only a limited set of parameters and thus limited freedom to fit them to input structures, it may be impossible to simultaneously reproduce two different physical quantities with the same parametrization of the potential.In such a situation, users must make a decision which physical properties are essential for their envisaged simulations and must choose or design a suitable parametrization of the potential matching their needs.Previously, a number of researchers published parametrizations [7][8][9][10][11][12][13][14] of this potential for elemental In, Ga, As, and the compound materials GaAs and InAs optimized for equilibrium bulk properties.Our detailed tests of these parameter sets with respect to the description of bulk elasticity and surface properties suggest that improvements are necessary ͑and possible͒ to accomplish quantitative investigations of the energetics of InAs/GaAs nanostructures.In this work, we present such an improved parametrization of the Abell-Tersoff potential for In, Ga, As, GaAs, and InAs.It simultaneously captures both bulk and surface properties with high overall accuracy.With this parametrization, it is now possible to perform reliable atomistic relaxations of InAs/GaAs nanostructures for quantitative studies of the total energies with dependable contributions from strained bulk and surface reconstructions.
Recently, Murdick et al. 15 presented a GaAs parametrization of an improved bond-order functional that additionally employs a penalty function to account for the electron counting rule for reconstructed surfaces.This very promising approach yields a good qualitative description of surface energies and defect energies and is able to describe intermediate bonding situations that might occur during, e.g., dynamic simulations of changes in the surface reconstruction.However, this capability is not necessarily needed in our simulations that are based on comparing the formation energy of nanostructures with given geometry.Furthermore, the potential of Murdick et al. 15 does not yet match the accuracy demands needed for our purposes: the relative deviation of the lattice constant ͑0.3%͒ and of the elastic constant c 44 ͑16%͒ of zinc blende GaAs may be too large for quantitative studies of lattice-mismatched nanostructures.To give a quantitative example, the elastic energy per volume ͓E el ͑2͒ ␣ 2 in Eq. ͑11͔͒ of GaAs that is biaxially strained in the ͕111͖ plane to the lattice constant of InAs is underestimated by about 13% using the results of Ref. 15.Our new parametrization ͑of the Abell-Tersoff functional͒ leads to a corresponding error of only 0.3%.Moreover, Ref. 15 provides a parametrization of GaAs only, whereas we require ͑and provide͒ potentials for both GaAs and InAs.
In Sec.II, we give the physical motivation and analytical form of the many-body potential as used in this work.Then, we describe the development of the parametrization in detail, provide the interaction parameters, and compare our parametrization to previous ones in terms of the properties of bulk phases and low-index surfaces.Section IV demonstrates the transferability of our potential to high-index GaAs and InAs surfaces and to Poisson ratios and elastic energies for various cases where biaxial strain is applied in a selected crystalline plane.We point out some of the limitations in describing defects and the diffusion of adatoms on surfaces.The remainder of this paper is then devoted to calculating the sizedependent formation energy of a representative QD and comparing the results to a previously employed hybrid approach that combined continuum-elasticity theory and DFT calculations ͑Refs.16 and 17͒.

II. ANALYTIC MANY-BODY POTENTIAL
The Abell-Tersoff potential was first introduced to model silicon. 18It was improved several times 5,19,20 and later on extended to the compounds SiC and SiGe ͑Ref.6͒ and many other systems.A similar functional form was introduced by Brenner 21 for the description of hydrocarbons.In the Abell-Tersoff potential, the total cohesive energy E coh of a configuration of atoms is given by a sum over all atoms i and their neighbors j, with the pairwise interatomic distance r ij = ͉r i − r j ͉.This expression is motivated by the suggestion of Abell 22 to describe the previously observed universal binding-energy curves for solid cohesion and chemisorption as a sum of next-neighbor pair interactions, which are functions of the local environment. 23The pairwise attractive and repulsive interactions are given by a Morse potential, The bond-order term B ij of the bond between atoms i and j depends on the neighborhood of both atoms: it weakens the attractive pair potential according to its other bonding partners.It is related to the energy difference of the bonding and the antibonding orbitals. 24The functional forms in the literature 6,9,11,12 slightly differ but can be brought to a common form by introducing the additional parameters n and m in the bond-order terms, The angular dependence of the bond order can be motivated from the second moment of the density of states, as obtained, e.g., from a tight-binding calculation, 24,25 and its particular form in the Abell-Tersoff potential corresponds to a bond. 26The angular function g ik ͑ ijk ͒ describes anisotropic interactions, characteristic for semiconductors with diamond or zinc blende structure, and is given by where ijk denotes the angle between r ij and r ik .The cutoff function f ij c ͑r ij ͒ limits the interaction range to atoms within a certain distance, In order to set up the broad list of reference data needed for the parameter optimization, we performed total-energy calculations in the framework of DFT.We use a plane-wave DFT code 42,43 together with norm-conserving ab initio pseudopotentials 44 and Monkhorst-Pack k-point meshes. 45lectron exchange and correlation is described in the localdensity approximation ͑LDA͒.For each physical property, we successfully performed convergence tests with respect to energy cutoff and number of k points.The cohesive energies E coh , lattice constants a 0 , and bulk moduli B were determined by fitting the equation of state to the Murnaghan equation of state.The surface calculations employed slabs and passivations with pseudohydrogen atoms, 46 as outlined in, e.g., Ref. 27.The surface energy ␥ per area A is given by the total energy of the slab E tot ͑after subtraction of its fixed lower part͒, the number of cations N III and anions N V , and the chemical potentials V and IIIV , 27 where the indices V and III denote As and Ga or In, respectively.In the following, we report the value of the surface energy ␥ for As =−E tot ͑As ␣ bulk phase͒, i.e., for a surface in equilibrium with bulk As, and for GaAs =−E tot ͑GaAs zinc blende bulk phase͒ and InAs =−E tot ͑InAs zinc blende bulk phase͒.The values for these total energies per atom of bulk structures are calculated with the same method as used for E tot in Eq. ͑8͒, i.e., by DFT-LDA calculations for the reference values of surface energies and by the Abell-Tersoff potentials if surface energies obtained with these potentials are quoted.
To avoid known problems of the LDA, e.g., the overestimation of cohesive energies, the reference data for the thermodynamic ground states of Ga, In, As, GaAs, and InAs were taken from experiment.The reference data for metastable bulk structures and for surfaces are not experimentally accessible and was taken from LDA calculations instead.Such a combination of theoretical and experimental data requires calibrating the structural trends of the DFT calculations with the available experimental measurements.To this end, we determined the ratios of experimental and LDA values for the lattice constant a 0 and for the cohesive energy E coh of the ground states.The LDA results without experimental counterparts were then scaled by these ratios in order to retain a consistent set of theoretical and experimental reference data.Instead of such a multiplication with the ratio of experimental and LDA values, the LDA data can be alternatively calibrated by a constant shift to the experimental value.While we find that some calibration procedure is required, the two alternatives gave comparable results, given that the input data can be matched by the parameter fitting only within a predetermined accuracy.

B. Parameter optimization
Our goal is to adapt the Abell-Tersoff potential to InAs/ GaAs nanostructures.In practice, this means that the param-  pose, we first selected a set of reference data representative for the bonding in the nanostructures and adjusted the parameters so as to minimize the error on this reference data set.Specifically, this data set consisted of selected bulk and surface properties ͑see Sec.III C͒.The error was calculated as the sum of weighted quadratic differences of the quantities, as obtained from the potential ͑with a certain set of param-eters͒ and the reference values.The error was minimized by applying the Levenberg-Marquardt algorithm 47 that smoothly varies from a steepest-descent minimization to an inverse-Hessian minimization as it approaches the minimum deviation value.Despite having many adjustable parameters, the flexibility of the Abell-Tersoff potential in describing different bonding situations is limited.This has several implications for parameter optimization: matching a reference property with its expression in terms of the potential constitutes a nonlinear equation in the potential parameters.Optimizing a set of potential parameters for a set of reference properties is thus equivalent to approximately solving a system of nonlinear equations.There is no general criterion if such a system is soluble ͑see, e.g., Ref. 47͒.In other words, from a mathematical point of view, it is not a priori known which subset of reference properties can be fitted with a requested accuracy.Hence, it is not necessarily the best choice to optimize all parameters for all reference data at the same time.There can be subsets of reference data that are conflicting in the sense that they cannot be simultaneously fitted with the requested accuracy.We overcame this obstacle by performing systematic optimization attempts that combine different subsets of parameters and different subsets of reference data.After identifying conflicting subsets of reference data and reducing the weight of one conflicting data, we repeated the search campaign to successively optimize the potential parameters.During this procedure, overfitting was avoided by checking against a test data set ͑properties of high-index surfaces, see Sec.IV A͒ not included in the optimization: parameter sets showing improvement only for the reference data, but not for test data, were discarded.Since the parameter fitting is a nonlinear minimization procedure, one might be worried about getting stuck in local minima.To cope with this risk, we sample the parameter space, using nearly 1000 different initial parameters sets to start the minimum search.These parameter sets were derived according to the procedure of Albe et al., 11 from previous parametrizations [7][8][9][10][11][12] or from successive fits of different subsets of parameters and reference data.

C. Numerical details
In this section, we describe how to calculate the physical quantities in our reference data from the Abell-Tersoff potential.The lattice constants a 0 and cohesive energies E coh were obtained by minimization of the cohesive energy E coh with respect to the lattice constant.The crystal lattices of ␣-Ga and ␣-As are determined by five ͑cf.Table II͒ and three ͑cf.Table IV͒ structural parameters, respectively.These structural degrees of freedom were individually determined by minimizing E coh while keeping the others fixed at their reference values.The bulk moduli B and elastic constants c 11 and c 12 were determined from numerical second derivatives of E coh : for B with respect to the volume and for c 11 and c 12 with respect to a variable in particular strain tensors. 48The elastic constant c 44 and the Kleinman parameter were calculated according to Nielsen et al. 28 The investigation of surfaces with a particular parametrization of the Abell-Tersoff potential requires to dovetail the lattice constant of the surface slabs with the value of the bulk lattice.For this purpose, we scale the DFT-LDA surface slabs to the value of the bulk lattice constant that we determined for a particular parametrization in the way described above.During the optimization of the parameters, this scaling was performed after every optimization step.We relaxed the atoms in the slab by a conjugate-gradient algorithm 47 until the maximum force in the system was below 1 meV/ Å.The surface energies for each parameter set were calculated by using Eq.͑8͒, inserting the cohesive energies as obtained with the same parameter set.
Furthermore, we defined two measures for the structural difference between the surfaces relaxed within DFT and with the Abell-Tersoff potential: first, with the Abell-Tersoff potential, we calculate the maximum atomic force in the input geometry ͑relaxed with DFT͒.This value, F 0 , would vanish if the geometry relaxed with DFT and with the Abell-Tersoff potential were identical.Note that a large absolute value of F 0 can also arise if, e.g., the bond length of a surface dimer is within the cutoff interval.Second, we calculate the average difference of all bond lengths ͗⌬r͘ in the surface unit cells after relaxation with either DFT ͑R ij DFT ͒ or the potential Note that the quantities F 0 and ͗⌬r͘ describe an extremal property and an average property, respectively, and therefore need not follow a simple proportionality.This is also observed in our calculations, indicating that the two measures are not redundant.

D. Parameters
The parameters for each pair of species as obtained from our optimization are given in Table I.In compliance with most previously published parametrizations, the parameters n ij and m ij were kept constant.The large value of D AsAs leads to an undesirable overestimation of the binding energy of the As dimer but was unavoidable to obtain proper surface energies: some of the considered surface reconstructions are terminated by As dimers, and these are only stable for large D AsAs .
Since In-Ga bonding is most likely not important for the InAs/GaAs nanostructures, such structures were not included in the reference data set.However, this set of parameters is sufficient to study ordering or segregation effects with a quasiequilibrium Monte Carlo simulation ͑see, e.g., Ref. 49͒ that is based on diffusion by cation exchange.For hightemperature molecular dynamics simulations of systems with both In and Ga atoms, the parameters for the possibly needed In-Ga interaction can be determined from an averaging scheme 6 that was shown to reproduce the cohesive energy and lattice constant of In x Ga 1−x As alloys. 12n the following, we compare the parametrization of the Abell-Tersoff potential for In, Ga, As, GaAs, and InAs developed in this work with previously published ones, which we chronologically name as T1, 7 T2, 9 T3, 8 T4, 10 T5, 11 T6, 12 T7, 13 and T8. 14 For each of these parametrizations, we use the identical bond-order functions ͓Eqs.͑4͒ and ͑5͔͒ as the authors by choosing n ij and m ij accordingly.For compound systems, we apply the same combinations of interaction parameter sets as in the original works.The parameter values for the interactions of In-In ͑T2͒, Ga-Ga ͑T1, T5͒, As-As ͑T1, T5͒, Ga-As ͑T1, T3, T5, T7, T8͒, and In-As ͑T2, T4, T6, T7, T8͒ are given in the corresponding publications.

E. Bulk structures
For each pairwise interaction, the cutoff radius truncates the range of the potential.The values were chosen similarly to previous parametrizations and limit the interaction to the nearest neighbors in In x Ga 1−x As zinc blende structures.Thus, for this choice of cutoff radii, the description of a zinc blende lattice with this potential is completely determined by the parameters of the Ga-As and In-As interactions.Nevertheless, the parameters for In-In, Ga-Ga, and As-As interactions are needed because E tot of the stable elemental crystal structures ͓i.e., ␣-Ga ͑A11͒, ␣-As ͑A7͒, and In͑bct͔͒ enters in the formula for the surface energy ͓see Eq. ͑8͔͒.A comparison of the results for the lattice parameters, cohesive energies, and bulk moduli of stable and metastable Ga, In, As, GaAs, and InAs bulk structures of experimental works, ͑scaled͒ DFT-LDA calculations, and calculations with our new parametrization and previous ones is given in Tables II-VI.For the internal structural degrees of freedom of Ga and As, we adopt the notation of Albe et al. 11 Some of the bulk moduli in Tables II-VI are given in brackets to indicate a technical artifact: equilibrium neighbor distances that are within the cutoff interval ͔ can lead to an unphysical influence of the curvature of the cutoff function ͓Eq.͑7͔͒ on the curvature at the minimum of the equation of state.This can cause incorrect results for quantities that are given by second derivatives of the total cohesive energy, such as the bulk modulus.
Table II shows that the parametrization developed in this work reproduces the lattice parameters and cohesive energies of all Ga structures with a few percent error.This implies that the Pauling relation between bond length and bond energy, as discussed in detail for T5 by Albe et al., 11 is fulfilled.In our optimization, the bulk moduli of the metastable structures turned out to form a conflicting subset with the GaAs surfaces.Therefore, we tolerated a larger error on these bulk moduli in favor of the ability to model the surfaces of nanostructures.
In Table III, we compile the investigated properties of In bulk structures.We find that they are reproduced with the same accuracy as the Ga bulk structures.The bulk moduli not explicitly included in the optimization procedure are described with a significantly higher accuracy, as compared to those of the Ga bulk structures.Note that the different quality in the description of In and Ga is due to the consideration of the InAs and GaAs surface properties in the optimization procedure.
The properties of As bulk structures ͑Table IV͒ are reproduced with slightly larger deviations, as compared to Ga and In.Here, one has to keep in mind that the As-As parameters were not only optimized with respect to As bulk structures but also to the properties of GaAs and InAs surfaces, where As surface dimers are a common structural element.Consequently, the As-As interaction is crucial for capturing the surface energies of these reconstructed surfaces.The dimerization of As leads to bond angles at the second-layer group-III atoms substantially different from those in the zinc blende structure.The angular terms of the III-As interactions thus penalize As-dimer formation and, therefore, a strong As-As interaction is required to overcompensate this penalty.Despite the resulting overestimation of the binding energy of the As dimer, the remaining As-As parameters still allow us to reproduce many bulk properties of arsenic.
The investigated structural and elastic properties of GaAs and InAs are compiled in Tables V and VI.Particularly, the properties of the zinc blende structures of these compound materials are reproduced with high accuracy by this potential.While most potentials give reasonable values for the lattice constant, cohesive energy, and bulk modulus, the present parametrization accurately reproduces all elastic constants, in particular, the shear modulus c 44 .The deviation of about 0.5 eV for the cohesive energy of the NaCl structure of GaAs was unavoidable in order to create a potential that simultaneously describes the reconstructed low-index GaAs surfaces.
These interactions are sufficient to describe heterogeneous systems of GaAs, InAs, and bulk interfaces formed by them, as well as In x Ga 1−x As alloys.Their applicability to subtle ordering effects in alloys still needs to be tested and is beyond the scope of this work.The description of point defects such as vacancies and interstitials is discussed in Sec.IV C.

F. Low-index compound surfaces
The surfaces of III-V semiconductors exhibit a noticeable number of different reconstructions.Most of them are actuated by the dimerization of surface As atoms.We include only a few typical reconstructions of low-index surfaces into the reference data set, namely, the ␣͑2 ϫ 4͒, ␣2͑2 ϫ 4͒, ␤͑2 ϫ 4͒, and ␤2͑2 ϫ 4͒ reconstructions of the ͑001͒ surface ͑cf.Fig. 1͒, as well as the ͑unreconstructed͒ ͑110͒ surface of GaAs and InAs ͑for the nomenclature, consult Ref. 36͒.For the growth of III-As nanostructures, moderately As-rich growth conditions are most important.As there is no indication that a coexistence of As-rich and cation-rich surface reconstructions has a strong impact on the growth process, we do not consider cation-rich surfaces, e.g., the ͑4 ϫ 2͒ reconstruction. 37The results obtained with our parametrization are compared to previous ones in Tables VII and VIII.The details of the considered properties and their calculation were given in Sec.III C.
This comparison clearly shows that the potential developed in this work achieves much higher overall accuracy in the description of the ͑110͒ and reconstructions of the ͑001͒ surfaces of both GaAs and InAs.All surface energies agree with the DFT-LDA reference values within about 10 meV/ Å 2 .Furthermore, the quantities F 0 and ͗⌬r͘ indicate that the agreement between the surface geometries as relaxed with DFT and with the analytic many-body potential is significantly improved, as compared to previous parametrizations.

IV. RELIABILITY TESTS
The so-called Stranski-Krastanov growth mode of semiconductor quantum dots ͑QDs͒ is driven by the balance between energy gain due to strain relief and energy cost due to the formation of side facets and edges. 17In this section, we will demonstrate that the above parametrization based on bulk properties and reconstructed low-index surfaces enables us to perform quantitative investigations of the formation energy and atomistic structure of lattice-mismatched InAs/ GaAs nanostructures such as QDs.In particular, it describes several high-index facets that are known to form side facets of QDs but also biaxially strained InAs films on GaAs with reasonable reliability.The limited applicability to bond breaking and making requires some caution when the potential is to be used for studies of defects and surface diffusion, as shown in Secs.IV C and IV D.

A. High-index compound surfaces
Recent high-resolution scanning-tunneling microscopy ͑STM͒ experiments revealed the atomic structure of InAs QDs on GaAs͑001͒ substrates 53,54 and on GaAs high-index substrates. 55,56Performing reliable studies of such systems with our potential calls for assessing the transferability from the low-index surfaces to which the potential was fitted to those GaAs and InAs high-index surfaces that are of particular importance as substrate and QD facets.We applied our potential to several high-index surfaces ͑cf.Fig. 2͒ with the same relaxation criteria as for the low-index surfaces and calculated the same quantities.A comparison of the results for GaAs and InAs high-index surfaces, as obtained from DFT-LDA calculations and from calculations with our parametrization and previous ones, is put together in Tables IX and X, respectively.
The error of the surface energy of all investigated highindex facets of about 10 meV/ Å 2 is very similar to the error for the low-index surfaces that were included in the optimization procedure.Note that every other parametrization gives deviations of more than a factor of 2 for at least one of the investigated low-index and high-index surfaces.Additionally, the maximum initial force on the DFT-relaxed highindex facets F 0 and the average bond-length deviation ͗⌬r͘ after relaxation with our potential are in the same range as for the low-index surfaces.This is a clear indication that our parametrization for low-index surfaces is transferable to the investigated high-index surfaces without loss of accuracy.The reason is that these high-index surfaces are reconstructed truncations of the zinc blende bulk dominated by As-dimer motifs, which also appear as the most prominent feature of the reconstructed low-index surfaces that were included in the fitting procedure.
Furthermore, we find that our parametrization reproduces the ͑4 ϫ 2͒ reconstruction of the GaAs͑001͒ and InAs͑001͒ surface within 5 and 20 meV/ Å 2 , respectively.However, the surface energies of reconstructions with an As bilayer ͓e.g., GaAs͑001͒c͑4 ϫ 4͒ ͑Ref.27͒ and In 2/3 Ga 1/3 As͑001͒␣2͑2 ϫ 3͒ ͑Ref.57͔͒ are overestimated by a factor of approximately 4. From the viewpoint of parameter optimization, these structures form a subset that is in conflict with the subset formed by the As bulk properties and the other surface reconstructions.Consequently, we did not include these two surface reconstructions as reference data in our optimization of the potential parameters.We note that none of the other investigated parametrizations gives acceptable results for both the ͑4 ϫ 2͒ reconstructions and the ones with an As bilayer.

B. Lattice-mismatched heterostructures
Biaxially strained bulk material is a good test case to check if this potential is able to capture the atomistic relaxation and elastic energy of epitaxial thin films, like the wet-ting layer that appears in Stranski-Krastanov growth of QDs.For the technologically relevant case of heteroepitaxial growth of InAs on GaAs substrates, the large lattice mismatch leads to significant elastic deformations in the InAs wetting layer.The term biaxial strain refers to a situation where a crystal is mechanically loaded in two directions, while it is free to relax in the third direction, which we identify by its Miller indices ͗hkl͘.The biaxial strain ␣, isotropic in each ͑hkl͒ crystal plane, causes an elastic relaxation ␤ along ͗hkl͘ related to ␣ by the biaxial Poisson ratio =−␤ / ␣ that mini-mizes the associated elastic energy E el .In the linear-response regime, the values of and E el of biaxially strained cubic crystal structures can be determined either analytically 58 from CET or numerically from atomistic simulations using strained simulation cells.Here, we employ both approaches to test the applicability of the potential developed in this work.
The analytic results were obtained by using the experimental elastic constants in the expressions derived from CET. 58 For the numeric results, we employed our potential in atomistic simulations of periodic supercells with two unit vectors of the strain plane and a set of strain values ␣ i ʦ ͓−0.07, 0.07͔.The minimization of the total energy E el ͑͒ for a certain biaxial strain ␣ i in the considered plane ͑hkl͒ with respect to the elastic response ␤ i perpendicular to this plane yields the biaxial Poisson ratio i =−␤ i / ␣ i .In each minimization, we relaxed the atoms until the relative change in total energy was less than 10 −3 .The elastic energy is obtained from the difference of the total energies of strained and unstrained supercells.The results of all investigated potential parametrizations show anharmonic strain dependencies and are well described by We can consistently assess the applicability of the potential developed in this work to the biaxial Poisson ratio and the elastic energy by comparing the harmonic terms ͑0͒ and E el

͑2͒
of the numerical results ͓Eqs.͑10͒ and ͑11͔͒ with the analytically obtained and E el .
For this comparison summarized in Tables XI and XII, we have exemplarily chosen biaxial strain in the ͑001͒ and ͑110͒ planes that correspond to extremal values of Poisson ratio and elastic energy, 58 and in the ͑137͒ and ͑3 7 15͒ planes as examples of two high-index planes.
Note that the formation of InAs quantum dots was reported only for GaAs substrate orientations that correspond to minimal or moderate elastic-energy densities of biaxially strained InAs, such as ͑001͒ and ͑113͒, 56 ͑114͒, 59 and ͑2 5 11͒, 55 but not for ͑110͒ 60 and ͑111͒ 61 substrates that would correspond to maxima in the elastic-energy density.In fact, STM experiments and ab initio calculations of heteroepitaxy of InAs on GaAs͑110͒ ͑Ref.62͒ suggested that the formation of dislocations is the dominant mechanism of strain relief.However, this competing mechanism of strain relief can be avoided by reducing the elastic-energy density, e.g., with a smaller value of the lattice mismatch ͓␣ in Eq. ͑11͔͒.Such a reduction with an additional InGaAs layer has recently proven to allow for QD growth even on GaAs͑110͒ substrates. 63he potential developed in this work reproduces the harmonic part of the biaxial Poisson ratio very well for all investigated strain planes in GaAs and InAs.The previously published parametrizations T3, T4, T6, and T7 give comparable results for the biaxial Poisson ratio and the elastic energy.Parametrizations T1 and T2 overestimate the biaxial Poisson ratio by a factor of 2 and show practically no dependence of the elastic energy on the applied biaxial strain.The weak dependence of ͑0͒ and E el ͑2͒ on the orientation of biaxial strain in the case of T8 is due to the overestimated value of c 44 in the original work. 50Previously published DFT calculations 58 showed that for InAs, the value of the anharmonic term E el ͑3͒ in Eq. ͑11͒ is of the order of −1 eV/ Å 3 for all Miller indices ͗hkl͘.This anharmonicity of the elasticenergy density is qualitatively reproduced by the InAs parametrizations T4, T6, and T7 and the one determined in this work.

C. Point defects
As a first example that reveals some of the limits of the Abell-Tersoff potential, we discuss neutral point defects in zinc blende GaAs with low formation energies ͑see, e.g., Ref. 64͒.These are vacancies ͑V As and V Ga ͒ and antisites with an As atom located at the lattice site of a Ga atom ͑As Ga ͒ or vice versa ͑Ga As ͒.Their formation energies as obtained from DFT calculations and with previous parametrizations of the Abell-Tersoff potential were given by Zollo et al. 64 Our calculations using the parametrization developed in this work and relaxations of simulation cells with 512 atoms are summarized in Table XIII.In the following, we report the point defect formation energies for As =−E tot ͑As ␣ bulk phase͒, i.e., for defects in equilibrium with bulk As.The formation energies of both vacancy defects, as obtained with the parametrization developed in this work, are improved over the previous ones, T1, T3, T5, and T8, and in reasonable agreement with the results from DFT calculations. 64The description of antisite defects is very bad, and for As Ga , we even obtain a negative formation energy.This is due to the occurrence of bonds among neighboring As atoms and the large value of D AsAs that was needed to stabilize surface reconstructions with As dimers.

D. Surface diffusion
As a second example that shows the limits of Abell-Tersoff potentials, we have chosen to investigate the energetics of bond breaking and making that is observable in potential-energy surfaces ͑PESs͒ of adsorbed atoms on surfaces.We exemplarily present the adsorption of a Ga atom on the ␤2 ͑2 ϫ 4͒ reconstruction of the GaAs͑001͒ surface, as obtained with the parametrization developed in this work.The results are compared to the DFT study by Kley et al. 40 and to a previous study by Salmi et al. 41 employing the T3 parametrization of the Abell-Tersoff potential.Additionally, we performed PES calculations with parametrization T5. 11 First, we relaxed a supercell containing 3 ϫ 5 units of reconstructed ͑2 ϫ 4͒ unit cells without adatom, keeping the lowest layers in the slab fixed.The resulting total energy is used to define the zero level of the adatom adsorption energy.Then, we determined the PES for Ga adsorption by positioning the Ga adatom on a grid of lateral positions to be held fixed and relaxing its height in the different lateral positions as well as the upper layers of the supercell until the absolute maximum force in the system was below 1 meV/ Å.The lateral spacing of the PES grid of ⌬x = ⌬y = a 0 / 25 results in more than 5000 points.This high resolution of the PES grid was chosen to avoid interpolation artifacts and to reveal localized artifacts that the cutoff function ͓Eq.͑7͔͒ may give rise to.The calculated differences of the total energies to that of the clean surface yields the PES of adatom binding energies shown in Fig. 3 for the parametrization determined in this work.The topology of the PES obtained with the parametrization determined in this work and with T5 are in qualitative agreement with the PES previously calculated using T3 ͑Ref.41͒: all PES show the four most stable adsorption sites, A1, A2, A3, and A4, and the overall corrugation of the PESs of about 2 eV is comparable to the corrugation of 1.75 eV from DFT calculations. 40Nevertheless, compared to the corresponding PES calculated with DFT, 40 all results with Abell-Tersoff potentials exhibit additional minima.We note that calculations of the surface diffusivity have been performed by computing the rate constants for hopping between adjacent minima within the approximation of transition state theory. 40ue to the additional minima and transition states encountered in the PES obtained with the Abell-Tersoff potentials, a direct comparison of hopping rates and energy barriers with those obtained in DFT is not informative.A meaningful comparison would require the calculation of the overall diffusivity, which is beyond the scope of the present work.The absolute depths of the realistic minima of the adsorption energy quantitatively depend on the parametrization, as shown in Table XIV.
The parametrization determined in this work somewhat underestimates the adsorption energy.The differences of the smallest and largest deviation relative to the DFT results are 0.71, 0.52, and 0.82 eV for T3, T5, and the parametrization developed in this work, respectively.The relative depth of the minima which is important for the adatom diffusion is not very well described by all these parametrizations.
To complete the picture of adsorption, we investigated the potential's ability to describe the effect of local changes on the hybridization.It was first observed by Kley et al. 40 that a Ga adatom has two stable adsorption geometries with either sp 2 -or sp 3 -like configuration at a surface As dimer of the GaAs͑001͒␤2͑2 ϫ 4͒ reconstruction.In similar calculations with the many-body potential developed in this work, a Ga  adatom in different heights above the As dimer and with different lengths of the latter exhibited only the sp 3 -like configuration and, thus, the analytic many-body potential does not capture this local change in the hybridization.The reason is that there is only one minimum in the angular function.We expect that introducing an additional minimum in the angular term of the Abell-Tersoff potential would allow the description of such multiple minima.

V. APPLICATION: QUANTUM DOT FORMATION ENERGY
Our parametrized Abell-Tersoff potential for the In-Ga-As system enables us to study the structure and energetics of realistic nanostructures.For a direct comparison of our atomistic method with previous works 16,17 that employed a hybrid method based on continuum-elasticity theory and DFT calculations, we determined the size-dependent formation energy of a freestanding InAs quantum dot ͑QD͒ on a GaAs͑001͒ substrate.An in-depth investigation of the relative stability of homogeneous films and differently shaped QDs will be presented elsewhere.

A. Details of the calculation
The following quantitative comparison with the hybrid approach is based on the simplest possible QD shape with exclusively ͕101͖ side facets.To investigate the variation of the QD formation energy, the calculations were performed for a series of QDs with different sizes.A meaningful comparison of both approaches requires a simple geometrical shape of the QD, i.e., the side facets must be full atomic layers.Adding layers increases each surface plane by a given orientation-dependent layer thickness.Therefore, in contrast to continuum approaches, an ideal isomorphic scaling is in many cases not achievable.In the special case of ͕101͖ facets, however, the QDs of different sizes are perfectly isomorphic.The number of In atoms in the QDs ranged from about 700 to about 20 000, corresponding to 10 and 31 atomic layers, respectively.
For the purpose of comparing both approaches, the QDs are assumed to sit on a homogeneous, 1.75 monolayer thick InAs wetting layer ͑WL͒ with the ␤2͑2 ϫ 4͒ surface reconstruction ͑Fig.4͒ that corresponds to moderately As-rich growth conditions.The detailed atomistic structure of the reconstructed WL and QD surfaces is based on DFT slab calculations of the corresponding free surfaces.Enforcing an overall pyramid shape of the QD creates low-coordinated atoms at the edges.While relaxing the atomic structure, these edge atoms are allowed to rearrange, but no attempt was made to enforce a particular atomic composition at the edges, e.g., by adding or removing low-coordinated edge atoms.Periodic boundary conditions in the ͑001͒ plane ensured that the lateral dimensions of the simulation cells were constant multiples of the GaAs lattice constant to model an ideal GaAs bulk underneath the InAs wetting layer.The lateral QD density that results from the area of our simulation cell corresponds to approximately 10 11 QD/ cm 2 , which is in line with typical experimental observations.We obviate finite-size effects by successfully performing convergence tests with respect to the vertical extension of the GaAs substrate.As a result, we find that the difference between the total energy of the chosen QD sizes and a homogeneous InAs film is converged to less than 1 meV per atom for substrates of more than approximately 10 nm thickness.The required vertical and lateral extension of the simulation cell results in systems of about 1.2ϫ 10 6 atoms.The individual structures of the QD series are relaxed with the manybody potential developed in Sec.II until the absolute value of the maximum force in the system was below 1 meV/ Å.In these slab calculations, the positions of the atoms in the lowest four atomic layers were kept fixed and all other atoms were allowed to move freely.

B. Results
The growth of self-assembled QDs is governed by the energetic balance of strain relief and the energetic costs of the formation of QD side facets and edges.Quantifying this interplay in terms of the QD formation energy is the key to understanding QD growth.For sufficiently large QDs, the contributions to the QD formation energy follow a scaling law 16 with respect to the QD volume V,

͑12͒
In our atomistic simulations, we numerically determine the formation energy E QD ͑V i ͒ for each relaxed structure i with QD of volume V i .This value is given as the difference between the total energy E tot QD of the slab with QD and the energy E tot sub of a slab with the wetting layer only.The different numbers of In and As atoms in the WL and the QD structures result in a dependency of the formation energy on the chemical potentials ͓cf.Eq. ͑8͔͒.Similar to the surface energy, the QD formation energy reads where ⌬E tot = E tot QD − E tot sub and ⌬N = N QD − N sub account for the difference in total energy and atom numbers, respectively.The comparison presented in the following is based on the formation energies for As =−E tot ͑As ␣ bulk phase͒, i.e., for a system in equilibrium with bulk As, and for InAs =−E tot ͑InAs zinc blende bulk phase͒.
For a direct comparison of our atomistic approach with the hybrid method of Pehlke et al., 16 we need to make con- tact between the number of atoms in our approach and the volume in the continuum representation of the QD in the hybrid approach.In the continuum model, the QD volume is unambiguously defined by the QD height h and scales with h 3 .For a consistent definition of the QD volume in our simulations, however, we need to additionally include atomistic contributions.In particular, we need to account for those atoms that complete the missing rows of the reconstructed WL surface, a contribution that scales with the QD base area, i.e., as h 2/3 .Therefore, a reliable leveling rule for our comparison of the atomistic and the continuum model is the QD volume but not the QD height.The QD volume in the atomistic simulations is directly related to the number of In and As atoms that form the QD.These can be determined by counting the atoms in a QD structure or by evaluating where n is the number of atomic layers that form the QD.The QD volume is then calculated as In our comparison, we have chosen the InAs lattice constant of the atomistic many-body potential ͑Table VI͒ to determine the volume per f.u.InAs V 0 = a 0 3 / 4. Calculating the formation energy ͓Eq.͑13͔͒ and the volume ͓Eq.͑15͔͒ of each QD structure in the isomorphic series ͑see Fig. 5͒ then enables us to determine the coefficients of the scaling law for the formation energy ͓Eq.͑12͔͒.The coefficients calculated from our atomistic simulations are in good agreement with the values previously obtained in the hybrid approach ͑see Table XV͒.The coefficient elastic accounts for the average strain energy per unit volume present in both the QD and the substrate after relaxation compared to unstrained InAs and GaAs bulk materials.When comparing the value of elastic , one has to keep in mind that the CET+ DFT value was obtained by using linear elasticity theory, while the atomistic approach includes nonlinear elastic effects.The deviations between our potential and the DFT results for the ͑unstrained͒ surface energies for the ͑001͒␤2͑2 ϫ 4͒ WL and the ͑101͒ QD side facets ͑Table VIII͒ mostly cancel each other.The deviation in the surface term is mainly due to the different dependency of the surface energy on biaxial strain.The edge contribution ⑀ edge accounts for undercoordinated atoms at the QD edges, but also for the fact that the unit mesh of the reconstruction, both on the side facets and in the wetting layer, is disturbed by the presence of edges.
The coefficients of the scaling law can be employed to quantify the range of validity of neglecting the edge contributions in the hybrid approach.In previous works employing the hybrid approach without the knowledge of the coefficient edge , this lower bound of the QD size could only be estimated.Our atomistically calculated coefficients, however, suggest that the relative weight of the edge contributions, i.e., the ratio edge V 1/3 / E QD ͑V͒, reaches a confidence limit of, e.g., 5% for QDs with more than approximately 2000 InAs f.u.For smaller QDs, the importance of the edge contributions increases significantly with decreasing QD volume.

VI. CONCLUSION
In this work, we provide a parametrization of the Abell-Tersoff potential that is particularly optimized for the description of InAs/GaAs nanostructures.It captures the lattice constants and cohesive energies of stable and several metastable bulk phases of In, Ga, As, GaAs, and InAs, as well as the elastic constants of GaAs and InAs with less than a few percent deviation.The energies of several reconstructed lowindex GaAs and InAs surfaces were simultaneously optimized with a remaining error of about 10 meV/ Å.We compared the atomic geometries for surface reconstructions obtained with DFT and with the many-body potential and optimized the potential parameters to yield improved geometries and energies compared to previous works.Our potential is transferable to the description of energies and structures of several reconstructed high-index GaAs and InAs surfaces.The harmonic parts of Poisson ratio and elastic energy of GaAs and InAs under biaxial strain in selected planes, as numerically calculated with our potential, are in TABLE XV.The coefficients of the scaling law for the QD formation energy ͓Eq.͑12͔͒, as obtained from atomistic simulations with the potential developed in this work, are in good agreement with the results of a previously employed hybrid approach ͑Ref.17, Fig. 7, filled square symbol͒.FIG. 5. ͑Color online͒ The QD formation energy, as obtained from our atomistic simulations ͑squares͒, follows the scaling law of ͓Eq.͑12͔͒ ͑black line͒ and is in good agreement with the hybrid approach ͑red line͒.The edge contributions are negligible for sufficiently large QD ͑dashed line͒.very good agreement with analytic results obtained from continuum-elasticity theory.The simultaneous description of bulk, surface, and elastic properties is a significant improvement to previously published parametrizations.The calculated defect formation energies and the potential-energy surface of a Ga adatom on a GaAs͑001͒␤2͑2 ϫ 4͒ surface show the typical limitations of Abell-Tersoff potentials.The calculated coefficients of the scaling law for the QD formation energy are in good agreement with a hybrid approach that employed strain fields from CET and surface energies from DFT.Moreover, the atomistic approach enables us to estimate the error in the hybrid approach due to the neglect of edge energies: for the particular atomic structure of the QD we investigated, the contribution from the formation of edges is less than 5% for QDs with more than approximately 2000 f.u.InAs.The parametrization developed in this work is suitable for quantitative studies of the energetic balance between strain relief and the formation of side facets and edges in defect-free InAs/GaAs heterostructures.In particular, it captures the strain tensors and total energies of overgrown and freestanding InAs QDs on GaAs substrates with good overall accuracy.

d
Value flawed due to influence of cutoff function.

FIG. 1 .
FIG. 1. ͑Color online͒ Reconstructions of the ͑001͒ surface of GaAs and InAs.Arsenic and cation atoms are shown in red ͑dark gray͒ and white, respectively.͑To guide the eyes, we show a 2 ϫ 2 repetition of the surface unit cells viewed along the crystallographic direction ͓001 ¯͔.͒

TABLE I .
Parameter sets for the different interactions.

TABLE II .
Equilibrium properties of Ga bulk phases as obtained from experiments and scaled DFT-LDA calculations, the potential developed in this work, and previous ones.
a Reference 11. b Reference 29.c Reference 30.d Value flawed due to influence of cutoff function.

TABLE III .
Equilibrium properties of In bulk phases as obtained from experiments and scaled DFT-LDA calculations, the potential developed in this work, and previous ones.
a Value flawed due to influence of cutoff function.b Reference 29.c Reference 30.

TABLE IV .
Equilibrium properties of As bulk phases as obtained from experiments and scaled DFT-LDA calculations, the potential developed in this work, and previous ones.

TABLE V .
Equilibrium properties of GaAs bulk phases as obtained from experiments and scaled DFT-LDA calculations, the potential developed in this work, and previous ones.
a Reference 31.b Values reported in original work overestimated ͑Ref.50͒.c References 32 and 33.d Reference 11. e Value flawed due to influence of cutoff function.

TABLE VI .
Equilibrium properties of InAs bulk phases as obtained from experiments and scaled DFT-LDA calculations, the potential developed in this work, and previous ones.
a Reference 34.b Values reported in original work overestimated ͑Ref.50͒.c Value obtained with phenomenological theory in Ref. 35. d Value flawed due to influence of cutoff function.

TABLE XI .
Poisson ratio for biaxial strain ͑0͒ determined analytically from CET ͑Ref.58͒ with the experimental elastic constants and numerically ͓Eq.͑10͔͒ with the parametrization determined in this work and with previous ones.

TABLE XII .
Elastic-energy density for biaxial strain E el ͑2͒ ͑eV/ Å 3 ͒ determined analytically from CET ͑Ref.58͒ with experimental elastic constants and numerically ͓Eq.͑11͔͒ with the parametrization determined in this work and with previous ones.

TABLE XIII .
Formation energies ͑eV͒ of vacancy and antisite point defects in zinc blende GaAs as obtained with DFT calculations ͑Ref.64͒ and the parametrization developed in this work.

TABLE XIV .
Adsorption energies ͑eV͒ of a Ga adatom at local minima A1, A2, A3, and A4 ͑Fig.3͒ on GaAs͑100͒␤2͑2 ϫ 4͒͒, as obtained with DFT, with the parametrization determined in this work and with previous ones.