text:
# Ferromagnetism in Mn doped GaAs due to substitutional-interstitial complexes

## I Introduction

The discussion theory ; cluster  of the physics that underlies room-temperature ferromagnetism in transition-metal doped semiconductors has largely focussed on substitutional geometries, e.g. MnGa site in GaAs. Indeed, there is a well-established tradition that 3$`d`$ impurities in III-V semiconductors are largely substitutional subst , while in Si they are mostly interstitials Ludwig . Modern first-principles total-energy calculations afford testing of this classic paradigm. Recent experiments rbs  find that Mn atoms occupy both substitutional as well as interstitial positions in GaAs. There have been suggestions from recent theoretical work erwin  suggested that primarily surface energetics will funnel Mn atoms in to interstitial sites from surface adatom positions. While MnGa behaves as a hole-producing acceptor, at the interstitial site, Mni behaves as an electron-producing donor. Since ferromagnetism is mediated by free-carriers, Mni could modify the magnetic properties from the case where only substitutional Mn sites were occupied.

Using density-functional theory as implemented within plane-wave pseudopotential total energy method, we consider here bulk and epitaxial growth conditions, investigating isolated defects (MnGa and Mni) and their complexes. We find that the Mn impurity in GaAs is stable in both substitutional and interstitial geometries depending on (a) the Fermi energy (which can be changed via external doping), (b) chemical potentials during growth and (c) bulk versus epitaxial growth conditions. The origin of these dependences is as follows: (a) As the formation energy of impurities that are neutral with respect to the lattice site they occupy (e.g. Mn$`{}_{Ga}{}^{}{}_{}{}^{0}`$) does not depend on the Fermi energy ($`ϵ_F`$), the formation energy of positively charged impurities (e.g. Mn$`{}_{}{}^{2+}{}_{i}{}^{}`$) decreases as $`ϵ_F`$ is shifted towards the VBM. Hence, the difference in the formation energies between Mn$`{}_{}{}^{2+}{}_{i}{}^{}`$ and Mn$`{}_{}{}^{0}{}_{Ga}{}^{}`$ decreases with $`p`$-type doping, resulting in increased solubility of interstitial Mn. (b) Substitution of Ga by Mn involves the removal of a Ga atom and the introduction of a Mn atom at the site vacated by Ga. Thus, substitution is generally enhanced under Ga-poor, Mn-rich growth conditions. On the other hand the formation energy of Mn at an interstitial site does not depend on the Ga chemical potential. Thus, one may stabilize substitutional (interstitial) doping using Ga-poor (Ga-rich) growth conditions. (c) Solid solubility can be controlled thermodynamically using epitaxial instead of bulk growth conditions wood . The absence of a substrate under bulk growth conditions allows the growing solid as well as its possible disproportionation products to attain their free-standing lattice geometry. This is the case when the growth takes place from the melt as in Bridgman growth. Then, if phase separation occurs, the precipitate will take up its most stable crystal structure, i.e. MnAs in the NiAs structure. In contrast, under thin-film epitaxial growth conditions (as in MBE, MOCVD) competing phases such as phase separated MnAs are forced to be coherent with the GaAs substrate. As zincblende MnAs strained on GaAs is less stable than the NiAs phase of MnAs, phase separation is more costly under coherent epitaxial conditions, and one expects wood  less phase separation, hence enhanced solubility. We find the following:

(i) Substitutional Mn has two stable charge states: neutral (Mn$`{}_{}{}^{0}{}_{Ga}{}^{}`$) and negatively charged (Mn$`{}_{}{}^{}{}_{Ga}{}^{}`$) charge state which have 1 and 0 holes, respectively. The calculated acceptor transition E(0/-) between these states occurs at Ev +0.13 eV, in good agreement with the experimental value of Ev+0.11 eV schneider . Here Ev corresponds to the valence band maximum of the host material.

(ii) The interstitial sites that Mn can occupy have either tetrahedral (coordinated to four As or four Ga atoms) or hexagonal symmetry. We find that Mn at the tetrahedral interstitial site coordinated by As is more stable than that coordinated by Ga, and exhibits a single charge state Mn2+ for all values of the Fermi energy. The (0/+) and (+/2+) donor transitions are found to lie inside the conduction band, so, isolated Mni produce electrons that will compensate the holes created by MnGa.

(iii) Under bulk growth conditions, the formation energy per Mn of substitutional Mn is $`\mathrm{\Delta }`$H(Mn$`{}_{Ga}{}^{}{}_{}{}^{0}`$)=0.91+$`\mu _{Ga}`$-$`\mu _{Mn}`$ eV=0.17 -$`\mu _{As}`$-$`\mu _{Mn}`$ eV, whereas interstitial Mn has $`\mathrm{\Delta }`$H (Mn$`{}_{i}{}^{}{}_{}{}^{2+}`$)=0.55-$`\mu _{Mn}`$+2$`ϵ_F`$ eV per Mn. Here $`ϵ_F`$ is the Fermi energy measured with respect to the valence band maximum of the host material, and $`\mu _{As}`$ and $`\mu _{Mn}`$ are the chemical potentials of As and Mn respectively. If we use maximally As-rich growth conditions ($`\mu _{As}`$=0 eV) and $`\mu _{Mn}`$=$`\mathrm{\Delta }`$H(MnAs), then we find $`\mathrm{\Delta }`$H(Mn$`{}_{}{}^{0}{}_{Ga}{}^{}`$)=0.91 eV. As one dopes the sample $`p`$-type and $`ϵ_F`$ approaches the VBM ($`ϵ_F`$=0), the energy difference between the formation energies of Mn$`{}_{}{}^{0}{}_{Ga}{}^{}`$ and Mn$`{}_{}{}^{2+}{}_{i}{}^{}`$ reduces to 0.38 eV. It could decrease even further if $`ϵ_F`$ penetrates the valence band with doping, or if the growth conditions are made less As-rich. The interstitial concentration is then expected to further increase. This is confirmed by recent experiments lpe  which use Ga-rich conditions for growth.

(iv) Under epitaxial growth conditions, the formation energies of both substitutional and interstitial Mn decrease by 0.74 eV/Mn, so their concentrations increase concomitantly leading to the possibilities of clusters. There is strong coulomb interactions between the oppositely charged constituents involving two substitutional (MnGa) and one interstitial (Mni); the cluster MnGa-Mni-MnGa is thus strongly stabilized and found to be more stable under $`p`$-type conditions than clusters involving three MnGa. Epitaxial growth conditions increases the solubility of such MnGa-Mni-MnGa clusters, with formation energy of -0.15 +$`ϵ_F`$ eV per cluster for the Q=+1 charge state under As-rich conditions and $`\mu _{Mn}`$=$`\mathrm{\Delta }`$H(MnAs).

(v) The presence of interstitial Mn in the MnGa-Mni-MnGa cluster provides a channel for the spins on the two substitutional Mn to align ferromagnetically even when there are no free carriers present in the cluster. We therefore conclude that ferromagnetism in GaAs:Mn can arise both from holes induced by isolated substitutional Mn atoms discussed previously (Dietl) as well as from charge compensated substitutional-interstitial clusters.

## II Method of Calculation

The formation energy for a defect comprising of atoms $`\alpha `$ in the charge state $`q`$ was computed using the density functional supercell method using the expression cuinse2

$`\mathrm{\Delta }H_f^{\alpha ,q}(ϵ_F,\mu )`$ $`=`$ $`E(\alpha )E(0)+{\displaystyle \underset{\alpha }{}}n_\alpha \mu _\alpha ^a+q(E_v+ϵ_F),`$ (1)
where $`E(\alpha )`$ and $`E(0)`$ are the total energies of a supercell with and without the defect $`\alpha `$ respectively. $`n_\alpha `$ denotes the number of atoms of defect $`\alpha `$ transferred in or out of the reservoir (equal to 1 for an atom removed and -1 for an atom added), while $`\mu _\alpha ^a`$ denotes their chemical potentials.

Total energies: The total energies of the charged supercells were computed by compensating any additional charge on the impurity atom by a uniform jellium background and have been corrected for interactions between charges in neighboring cells using the Makov and Payne correction payne . For isolated defects we used both the monopole as well as quadrupole corrections, while for composite defects we have added only the monopole correction to the total energy assuming all the charge to be localized at a single point. We use the static dielectric constant of GaAs (12.4) madelung . The quadrupole moment of the isolated defects was calculated as the difference between the moments of the supercell with the charged defect and that with the neutral defect.

Transition energies: The defect transition energy $`ϵ(q,q^{})`$ is the value of the Fermi energy $`ϵ_F`$ at which $`\mathrm{\Delta }H^{\alpha ,q}(ϵ_F)`$=$`\mathrm{\Delta }H^{\alpha ,q^{}}(ϵ_F)`$. The zero of the Fermi energy is chosen as the valence band maximum $`E_v`$ of the pure host at the $`\mathrm{\Gamma }`$ point.

Chemical potential limits: As the reservoir supplying the atoms could be elemental solids, or compounds formed from the elements, we express $`\mu _\alpha ^a`$ as the sum of the energy of the element in its most stable structure $`\mu _\alpha ^s`$, and an additional energy $`\mu _\alpha `$ i.e $`\mu _\alpha ^a`$ = $`\mu _\alpha ^s`$ \+ $`\mu _\alpha `$. The required ranges of $`\mu _\alpha `$ are determined by $`\mu _{Ga}`$ $``$ 0; $`\mu _{Mn}`$ $``$ 0; $`\mu _{As}`$ $``$ 0 (no precipitation of solid elements) and by the formation energies of GaAs, and MnAs. The allowed values of chemical potential are such that GaAs is stable, i.e. $`\mu _{Ga}+\mu _{As}=\mathrm{\Delta }H_f(\text{GaAs})`$, the latter being the formation energy of zinc-blende GaAs. Further, as Mn should not precipitate as MnAs, we restrict $`\mu _{Mn}+\mu _{As}<\mathrm{\Delta }H_f(\text{MnAs})`$, the formation energy of MnAs in its most stable (NiAs) structure. For epitaxial growth conditions, the formation energy of zinc-blende MnAs lattice-matched to GaAs is considered. In this case we calculate the epitaxial formation energy, $`\mathrm{\Delta }H_f(\text{MnAs})_{epi}`$, forcing the in-plane lattice constant of MnAs to equal that of GaAs, while the out-of-plane lattice constant, c, is allowed to vary. For coherent epitaxial growth the condition that MnAs should not form during incorporation of Mn in GaAs becomes $`\mu _{Mn}+\mu _{As}<\mathrm{\Delta }H_f(\text{MnAs})_{epi}`$.

The energies E($`\alpha `$), E(0), $`\mathrm{\Delta }H_f(\text{GaAs})`$, $`\mathrm{\Delta }H_f(\text{MnAs})`$, $`\mathrm{\Delta }H_f(\text{MnAs})_{epi}`$ and $`\mu _\alpha `$ are calculated within the density functional formalism, through the momentum-space pseudopotential total energy representation ihm , using ultrasoft pseudopotentials ultrasoft . The GGA-PW91 version of the exchange-correlation functional pw91  and no correction for the band gap underestimation was made. The calculations were performed over a Monkhorst-Pack 4x4x4 k-point grid for 64 compare  and 216 atom supercells of GaAs using VASP vasp . Changing the k-point mesh from 2x2x2 to 4x4x4 changed the formation energies by $``$ 20 meV. Larger 256 atom supercells with 1x1x2 k-points were used for the calculations with clusters to ensure a larger separation between clusters. We used a plane wave cutoff of 227.2 eV for these calculations. Increasing the cutoff to 300 eV, changed the formation energies by $``$ 10 meV. As the lattice constant of the supercell was kept fixed at the GGA optimised value for GaAs of $`a`$=5.738 $`\AA `$ footnote1 , the internal coordinates were optimised. Our calculated (experimental) formation energies are $`\mathrm{\Delta }H_f(GaAs)`$ = -0.74 (-0.74), $`\mathrm{\Delta }H_f(MnAs)`$ = -0.74 (-0.61) eV and $`\mathrm{\Delta }H_f(MnAs)_{epi}`$ $``$ 0 eV. For elemental Mn we assume the nonmagnetic fcc structure mn\_reference , while for elemental Ga, we assume the base-centered orthorhombic structure.

The charge corrected payne  MnGa (0/-) transition as well as the difference in formation energies between Mn$`{}_{}{}^{0}{}_{Ga}{}^{}`$ and Mn$`{}_{}{}^{2+}{}_{i}{}^{}`$ are given in Table I for supercell sizes of 64 and 216 atoms. We see that changing the supercell size from 64 to 216 atoms lowers the acceptor energy by 30-50 meV and stabilizes Mn$`{}_{}{}^{0}{}_{Ga}{}^{}`$ over Mn$`{}_{}{}^{2+}{}_{i}{}^{}`$ by 50-150 meV. The charge correction increases the acceptor energy by 60-90 meV and stabilizes Mn$`{}_{}{}^{0}{}_{Ga}{}^{}`$ over Mn$`{}_{}{}^{2+}{}_{i}{}^{}`$ by 250-350 meV.

## III Results

### III.1 Isolated substitutional Mn on the Ga site of GaAs

Fig. 1 describes the formation energy $`\mathrm{\Delta }`$H(Mn$`{}_{}{}^{0}{}_{Ga}{}^{}`$) of neutral substitutional Mn in GaAs as a function of the chemical potentials $`\mu _{As}`$ and $`\mu _{Mn}`$. The shaded areas denote chemical potentials that produce unwanted products: (i) When $`\mu _{As}`$ becomes greater than zero (the cohesive energy of solid As) we have precipitation of elemental As as shown on the left hand side of Fig. 1. (ii) In the opposite limit, when $`\mu _{As}`$ takes more negative values than the formation energy $`\mathrm{\Delta }`$H(GaAs), we have maximally As-poor conditions and the host itself becomes unstable, as shown on the right hand side of Fig. 1. (iii) The diagonal lines in the main body of Fig. 1 denote different values of $`\mu _{Mn}`$. When the chemical potential of Mn becomes greater than zero (the cohesive energy of solid Mn), metallic Mn will precipitate as shown in the bottom right corner of Fig. 1. Conversely, (iv) when $`\mu _{Mn}`$ becomes equal or larger than $`\mathrm{\Delta }`$H(MnAs)-$`\mu _{As}`$, we will precipitate a secondary phase of MnAs. Clearly, since $`\mu _{Ga}`$+$`\mu _{As}`$=$`\mathrm{\Delta }`$H(GaAs) and $`\mu _{Mn}`$+$`\mu _{As}<`$ $`\mathrm{\Delta }`$H(MnAs), one can keep the latter inequality even for moderately negative values of $`\mu _{Mn}`$, provided that $`\mu _{As}`$ is adjusted. The lines in Fig. 1 show that the lowest $`\mathrm{\Delta }`$H(Mn$`{}_{}{}^{0}{}_{Ga}{}^{}`$) value is 0.91 eV (circle at bottom left corner). This can be attained at $`\mu _{As}`$=0 (maximally As-rich); $`\mu _{Mn}`$=$`\mathrm{\Delta }`$H(MnAs). Alternatively, the same solubility can be attained for less rich-As conditions, but richer Mn conditions, e.g. for $`\mu _{As}`$=-0.5 eV and $`\mu _{Mn}`$=-0.24 eV.

Having described in Fig. 1 the stability of the neutral substitutional, we next describe in Fig. 2 the stability of the charged substitutionals. Here we chose the chemical potentials $`\mu _{As}`$=0, $`\mu _{Mn}`$=$`\mathrm{\Delta }`$H(MnAs) (denoted by the circle in Fig. 1) and vary the Fermi energy. We see that for $`p`$-type conditions, the lowest energy charge state is Mn$`{}_{}{}^{0}{}_{Ga}{}^{}`$, whereas for higher Fermi energy the stablest charge state is Mn$`{}_{}{}^{}{}_{Ga}{}^{}`$. Table I gives the (0/-) acceptor transition energy calculated with various supercell sizes with and without charge correction. The most converged (0/-) transition energy calculated for the 216 atom cell and corrected for charge interactions is Ev+0.13 eV, in good agreement with the measured value of Ev+0.11 eV schneider . Fig. 2 shows that under epitaxial conditions (right y-axis), the formation energy of Mn$`{}_{}{}^{0}{}_{Ga}{}^{}`$ is lowered by 0.74 eV.

We next describe the electronic structure of MnGa. In Figs. 3(a) and (b) we show the Mn $`d`$ projected partial density of states (PDOS) for two charge states of substitutional Mn. The main features can be understood as arising from the hybridization between the anion dangling bonds generated by a Ga vacancy and the $`d`$ levels on the Mn ion placed at the vacant site subst . The Mn $`d`$ ion levels are split by the tetrahedral crystal field into $`t_2(d)`$ and $`e(d)`$. Exchange interactions further split these levels into spin-up ($``$) and spin-down ($``$) levels. The $`t_2(d)`$ levels on the Mn atom hybridize with the levels with the same symmetry on the As dangling bonds, while the $`e(d)`$ levels have no other states available for significant coupling subst . Because the location of the Mn ion $`d`$ levels is below the dangling bond levels, after hybridization, the deeper bonding $`t_2`$ states have dominantly Mn $`d`$ character (refered to as CFR: ”crystal field resonance”), while the higher antibonding $`t_2`$ states have dominantly As $`p`$ character (refered to as DBH : ”dangling bond hybrid”). These interactions lead to the energy level diagram depicted schematically on the left-hand-side of Fig. 4 showing a fully occupied, Mn-localized up-spin CFR of $`t_2`$ and $`e`$ symmetries. At a higher energy we have the up and down-spin DBH states with $`t_2`$ symmetry. Because of the location of the Ga vacancy states $`t_2(p)`$ between the exchange split $`t_2(d)`$ states on the Mn, a negative exchange splitting is induced as a result of hybridization on the DBH states our\_paper  with $`t_{DBH}^{}`$ below $`t_{DBH}^{}`$. As a result, the neutral substitutional defect Mn$`{}_{Ga}{}^{}{}_{}{}^{0}`$ has the electron configuration \[$`t_{}^3e_{}^2`$\]CFR ($`t_{}^3t_{}^2`$)DBH, with a total magnetic moment $`\mu `$=4 $`\mu _B`$, and a hole in the $`t_{DBH}^{}`$ orbital. This configuration corresponds to the multiplet 5T2 as observed in electron paramagnetic resonance (EPR) experiments mn\_acc . The partial occupancy of the negative exchange-split DBH states stabilizes the ferromagnetic state over the antiferromagnetic state our\_paper .

### III.2 Isolated interstitial Mn

Mn interstitial can occupy a site with tetrahedral symmetry (coordinated by four As or four Ga atoms) or a site with hexagonal symmetry. We have calculated the total energies of Mn at these positions in a 64 atom cell of GaAs, and the results for the tetrahedral interstitial sites are given in Table II. The tetrahedral interstitial Mni(As) coordinated by four As atoms is more stable than the one coordinated by four Ga atoms, with the difference being 0.16, 0.31 and 0.31 eV for charge states $`q`$=1,2 and 3. In contrast, the hexagonal interstitial has 0.62 eV higher total energy than the most stable Mn$`{}_{i}{}^{}{}_{}{}^{2+}`$(As). Experimentally, the presence of interstitial Mn was detected by an analysis of the EPR spectrum epr  as well as by Rutherford back scattering rbs . The distinction between the two types of $`T_d`$ interstitial sites (Mn-next to As vs Mn-next to Ga) is difficult to determine experimentally and involved an analysis of the experimentally measured contact interaction in terms of the covalency of the Mn-X bond. This analysis suggested that Mni(Ga) was more stable, while our total energy calculations suggest that Mni(As) is more stable.

The formation energy of various charge states of interstitial Mn is shown in Fig. 2 for $`\mu _{As}`$=0 and $`\mu _{Mn}`$=$`\mathrm{\Delta }`$H(MnAs). We see that the stable charge state is Mn$`{}_{}{}^{2+}{}_{i}{}^{}`$ for the full range of Fermi level, with maximum stability at $`ϵ_F`$=0. To compare the relative stability of Mn$`{}_{}{}^{2+}{}_{i}{}^{}`$ at $`ϵ_F`$=0 with substitutional Mn$`{}_{}{}^{0}{}_{Ga}{}^{}`$, we show in the upper scale of Fig. 1 the difference $`\mathrm{\Delta }`$H(Mn$`{}_{}{}^{2+}{}_{i}{}^{}`$) - $`\mathrm{\Delta }`$H(Mn$`{}_{}{}^{0}{}_{Ga}{}^{}`$) between the formation energies of interstitial and substitutional Mn. We see that substitutional Ga is stabler on the left hand side of the figure, i.e sufficiently As-rich, whereas interstitial Mn is stabler at the right hand side of the figure, i.e. sufficiently As-poor. The energy difference is

$$\mathrm{\Delta }H(Mn_i^{2+})\mathrm{\Delta }H(Mn_{Ga}^0)=0.38+\mu _{As}+2ϵ_F$$

For $`\mu _{As}`$=0, the substitutional Mn are stabler by 0.38 eV, while for moderately As-rich conditions, say $`\mu _{As}`$ = -0.4 eV, both defects have comparable formation energies.

These results are in agreement with recent experiments using liquid phase epitaxy lpe  to introduce Mn in GaAs. Experimentally a decrease in hole concentration is found as the Mn concentration is increased. Under the Ga-rich growth conditions used, As antisites are not expected to be the dominant source of the observed compensation. Hence the major source of compensation is believed to come from Mni as expected for Ga-rich conditions from Fig. 1.

We next examine the electronic structure of Mn at a tetrahedral interstitial site. When Mn occupies a tetrahedral interstitial position, five of the seven electrons occupy the $`t_{}e_{}`$ CFR levels, with the remaining two going into the down-spin $`t_{}`$ levels. This is evident from the PDOS for the doubly-ionized Mn$`{}_{}{}^{2+}{}_{i}{}^{}`$ shown in Fig. 3(c), where Mn$`{}_{}{}^{2+}{}_{i}{}^{}`$ is found to have the configuration \[$`t_{}^3e_{}^2`$\]CFR with a magnetic moment of $`\mu `$=5$`\mu _B`$. The central panel of Fig. 4(a) shows schematically the levels of Mni. As the (0/+) and (+/2+) transitions are calculated to lie inside the GaAs conduction band (Fig. 2), we conclude that Mni produce free electrons in GaAs.

### III.3 Clusters of substitutional Mn

Having dealt with the isolated limit, we investigated whether Mn atoms show a tendency to cluster. Recent experiments cocluster  on dilute magnetic semiconductors have found a strong tendency of the doped transition metal atoms to cluster and there has been some theoretical work cluster  to support such observations. We consider As-centred clusters \[(As)MnnGa4-n\] with $`n`$=0,1,2,3 and 4.

Fig. 5 shows the formation energy of clusters made of three substitutional Mn atoms (S-S-S) at lattice locations (0,0,0), (a/2,a/2,0) and (0,a/2,a/2) in the 64 atom supercell. This corresponds to the $`n`$=3 cluster. Here $`a`$ is the cubic lattice constant of GaAs. We see that the neutral cluster (3 MnGa)0 having three holes is stable under $`p`$-type conditions, whereas the charged cluster (3MnGa)- with 2 holes is more stable above $`ϵ_F`$ = 0.15 eV. The energies of the complex with 3(MnGa)0 is 2.2 + 3($`\mu _{Ga}`$ \- $`\mu _{Mn}`$) eV while that of three noninteracting MnGa in their lowest energy charge state is 2.71 + 3($`\mu _{Ga}`$ \- $`\mu _{Mn}`$) eV. For epitaxial conditions $`\mathrm{\Delta }`$Hepi=0.02 +3($`\mu _{Ga}`$ \- $`\mu _{Mn}`$) eV/cluster. Thus, as the formation energy is very low, the tendency for the Mn atoms to cluster is strongly enhanced under epitaxial growth conditions.

In order to obtain a measure of the tendency to cluster, we calculate the clustering energy. The ”clustering energy” $`\delta (n)`$ is defined as the energy difference between n substitutional Mn atoms surrounding an As site \[(As)MnnGa4-n; $`0n4`$\] and n isolated well-separated constituents. Thus, $`\delta E(n)=[E(n)E(0)]n[E(1)E(0)]`$, where E(n) is the total energy of the supercell with As-centered clusters of n Mn atoms. We find that $`\delta E`$(n)= -228, -482 and -794 meV per cluster of n=2,3 and 4 Mn atoms for a 64 atom supercell. The clustering energy changed to -519, -1069 meV for clusters involving 2 and 3 Mn in a 256 atom supercell. These results indicate a strong tendency for the neutral substitutional Mn atoms to form clusters.

### III.4 Neutral complexes of MnGa and Mni

We considered the defect complex formed between MnGa (S) and interstitial Mni (I) denoted as (S-I-S)Q, where $`Q`$ is the total charge of the complex. The geometry for the complex had MnGa at (0,0,0) and (a/2,a/2,0) and Mni at (a/2,0,0) in the 64 atom supercell of GaAs. We see in Fig. 5 that S-I-S exists in two charge states: When the Fermi energy is below Ev+0.1 eV we have the stable structure is (S-I-S)1+, whereas when $`ϵ_F`$ is above it, the stable structure is the neutral (S-I-S)0. Thus, the donor transition for the cluster is at Ev+0.1 eV. Fig. 5 also shows that for Fermi levels below Ev+0.22 eV, the S-I-S complex is stabler than the S-S-S complex. As for the interaction energy of the components of the complex: the formation energy of the non-interacting neutral ($`Q=0`$) components of the complex is 2E(Mn$`{}_{Ga}{}^{}{}_{}{}^{0}`$)+E(Mn$`{}_{}{}^{0}{}_{i}{}^{}`$)=4.31+2$`\mu _{Ga}`$ \- 3$`\mu _{Mn}`$ eV per 3 impurities, while the formation energy of the interacting neutral complex is 1.41+2$`\mu _{Ga}`$-3$`\mu _{Mn}`$ eV. This represents a $``$ 2.9 eV per 3 impurities stabilization over the non-interacting, neutral defects. The energy of the neutral complex measured with respect to the stablest lattice site occupied by isolated Mn under particular experimental conditions is found to be -574 meV for $`\mu _{As}`$=0 eV, $`\mu _{Mn}`$=$`\mathrm{\Delta }`$H(MnAs) and $`ϵ_F`$=0 eV. Hence this complex is strongly stabilized.

The reasons for the stability of the (S-I-S)0 complex can be appreciated from Fig. 4(a). Upon bringing together 2Mn$`{}_{}{}^{0}{}_{Ga}{}^{}`$ with Mn$`{}_{}{}^{0}{}_{i}{}^{}`$, one electron drops from the higher energy $`t_{}^{CFR}`$ level of Mni to the lower energy $`t_{}^{DBH}`$ level of each substitutional site, resulting in $`[t_{}^3e_{}^2]_{CFR}(t_{}^3t_{}^3)_{DBH}`$ configuration at each MnGa site (Fig. 4(b)) which corresponds to Mn$`{}_{}{}^{}{}_{Ga}{}^{}`$. These conclusions are evident from our calculated DOS of the S-I-S complex, projected on the I and S sites shown in Fig. 6. We find that for both $`Q`$=0 (Fig. 6(a)) and $`Q`$=2 (Fig. 6(b)) the I site has the configuration $`[t_{}^3e_{}^2]_{CFR}`$ or ”$`d_{}^5`$”. This substitutional-to-interstitial charge transfer lowers the energy of the complex by twice the separation between $`t_{}^{CFR}`$ level of Mni and $`t_{}^{DBH}`$ level of MnGa. Furthermore, it creates a favorable Coulomb attraction between the components S--I2+-S- of the complex. This energetically favorable substitutional-interstitial association reaction then eliminates the holes that were present in isolated substitutional MnGa and could explain the puzzling observation hole\_conc  of the existence of a far lower concentration of holes than Mn in GaAs. Alternate explanations such as the presence of As antisites hole\_mech  as well as the presence of Mn atoms connected to six As atoms (as in the NiAs structure) have been offered. However, samples have been prepared where the concentration of As antisites is too low to explain the observed compensation of holes. Further, six-fold coordinated Mn atoms have not been observed in Mn doped GaAs samples exafs .

### III.5 Ferromagnetism of the $`(SIS)^0`$ complex

The neutral complex has two Mn$`{}_{}{}^{}{}_{Ga}{}^{}`$ and one intervening ”$`d_{}^5`$” interstitial. We find that a ferromagnetic arrangement between MnGa is favored in the complex ($`Mn_{Ga}^{}Mn_i^{2+}Mn_{Ga}^{}`$)0 by 176 meV. In contrast, our calculations for two Mn$`{}_{}{}^{}{}_{Ga}{}^{}`$ atoms without the intervening interstitial atom finds that an antiferromagnetic arrangement of spins on the substitutional Mn atoms is favored by 108 meV. Thus Mni is responsible for mediating a ferromagnetic interaction between Mn$`{}_{}{}^{0}{}_{Ga}{}^{}`$.

How does the presence of the interstitial Mn mediate the alignment of spins on the substitutional Mn? There are three possible arrangements for the spins on the Mn atoms making up the neutral complex \- (SIS)0, (SIS)0 and (SIS)0. From our total energy calculations we find that the energy for the configurations (SIS)0 and (SIS)0 are higher by 563 meV and 176 meV, respectively, than the energy, E0, of the ground state (SIS)0. (The energies changed marginally to 602 meV and 192 meV respectively when we increased the supercell size to 256 atoms.) The stabilization of the SIS magnetic arrangement can be understood using simple arguments: In the configuration (SIS)0, as one spin channel is completely filled, there is no channel of hopping available for the electrons to delocalize and lower their energy. Thus, this is a high energy spin configuration with energy E0+563 meV. In contrast, in the configuration (SIS)0, two channels of hopping are present; the first between the electrons on S and S, and the second between those on S and I. This configuration is found to have the energy, E0+176 meV. Likewise, the configuration (SIS) which has energy E0 has two channels of hopping present between S and I. The dominant factor in determining the configuration which has the lowest energy are the hopping matrix elements - VS,I between S and I and VS,S between the two S. To a first approximation, these hopping matrix elements are determined by the separation between the atoms involved. As the distance between the two substitutional Mn atoms is $`\sqrt{2}`$ times the distance between S and I, the effective hopping matrix element between S is smaller. Hence the presence of an intervening Mni provides a channel for the ferromagnetic arrangement of spins between two MnGa even in the neutral charge-compensated complex. In contrast, the presence of a closed shell donor such as AsGa between two MnGa gives rise to an antiferromagnetic (or weakly ferromagnetic) interaction between MnGa.

How does the presence of the interstitial affect long-range ferromagnetism? In order to investigate this we introduced a hole-producing, isolated substitutional Mn atom at different lattice locations, indicated in Fig. 7, and investigated whether the spin on this isolated substitutional Mn atom prefers to align parallel or antiparallel with respect to the spins on the substitutional Mn atoms within the S-I-S cluster. We find that the substitutional Mn likes to align ferromagnetically with the MnGa of the cluster by 147, 214 and 81 meV, respectively for the positions 1, 2 and 3 (see Fig. 7). Hence the presence of the interstitial Mn forces a hole which is located $``$ 12 $`\AA `$ from the S-I-S cluster to align ferromagnetically and therefore contributes to the long-ranged ferromagnetism observed in these systems.

As discussed earlier, the basic electronic structure of substitutional Mn in GaAs can be understood as arising from the hopping interaction between the Mn $`d`$ states and the As $`p`$ dangling bond states. Therefore, the coupling between two Mn atoms is through the As $`p`$ states. It is strongest along the directions in which the ”$`pd`$” coupling of the Mn with the As states is the largest and decreases with distance along that direction. When there is a hole, the antibonding $`t_2^{}`$ orbitals are partially occupied, and there is ferromagnetism. Hence in Fig. 7 spins on the Mn atoms at sites 2 and 3 prefer to align ferromagnetically with the spins on S when there is partial compensation. This could happen either because the S-I-S cluster is totally compensated, and a hole exists on the Mn at sites 2/3 or the S-I-S cluster is partially compensated. The mechanism stabilizing the ferromagnetic coupling between S and site 1 is the same as what we discussed for the (S-I-S)0 complex earlier and exists even when there is total compensation.

### III.6 Charged MnGa-Mni-MnGa complexes:

While the neutral complex has no holes, the $`Q`$=+1 and +2 complexes have 1 and 2 holes respectively (Fig. 4(c)) and a net magnetic moment of 4 and 3 $`\mu _B`$ respectively. For $`Q`$=+2, MnGa adopts the configuration Mn$`{}_{}{}^{0}{}_{Ga}{}^{}`$ (Fig. 4(c)). We find that (MnGa-Mni-MnGa)2+ prefers the ferromagnetic arrangement of spins on MnGa by 286 meV, similar to the ferromagnetic preference ($``$ 305 meV) of Mn$`{}_{}{}^{0}{}_{Ga}{}^{}`$-Mn$`{}_{}{}^{0}{}_{Ga}{}^{}`$ pair without an intervening Mni. These results suggest the surprising fact that the spins on MnGa align ferromagnetically in the charged complexes, almost as if Mni did not exist. The number of holes in the cluster is the same as the number in the pair, though the number of Mn atoms are different. This is in agreement with the experimental observation mag\_deficit  where above a critical concentration of Mn, both the number of holes as well as the ferromagnetic transition temperature remain constant , while the magnetic moment per Mn atom decreases mag\_deficit . The magnetic moments that we obtain for the $`Q`$=+1 and +2 charge states translate into average moments of 1.33 $`\mu _B`$ and 1 $`\mu _B`$ per Mn, while the uncompensated pair of MnGa have a magnetic moment of 4 $`\mu _B`$. In this regime where the $`T_c`$ is found to saturate, the average magnetic moment per Mn is found to vary from $``$ 3 $`\mu _B`$ at a Mn concentration of 5.54 $`\%`$ to 1.74 at 8.3 $`\%`$.

Recent experiments waluck2  find that the $`T_c`$ of the as-grown samples increased after annealing. This was interpreted as the migration of FM-reducing interstitial Mn to FM-enhancing substitutional positions. We investigated which clusters could break by annealing and promote ferromagnetism. As the S-I-S complexes are rather strongly bound with respect to their constituents, we investigated instead complexes S-I, which are bound weakly ($``$ -196 meV in the +1 charge state for $`\mu _{As}`$, $`\mu _{Mn}`$=$`\mathrm{\Delta }`$H(MnAs) and $`ϵ_F`$=0 eV). We find antiferromagnetic spin arrangement in all $`Q`$=0, +1, +2 and +3 charge states considered. Thus, when these weakly bound S-I clusters are broken, depending on the charge state, there could be an increase in the number of holes and consequently the ferromagnetic transition temperature. On the other hand, S-I-S clusters appear to be stable and hence do not disintegrate under annealing.

## IV Summary

Under As-rich conditions, Mn prefers to substitute the Ga site. As the growth conditions become less As-rich, or as extrinsic doping pushes $`ϵ_F`$ towards and even below the VBM, the formation energy of interstitial Mn becomes competitive with that of substitutional Mn. Under coherent epitaxial growth conditions, when MnAs precipitates are forced to be coherent with the zinc-blende lattice, the formation energy of both substitutional and interstitial decrease. At this point, the solubility is large enough to form clusers. We find that S-I-S clusters are more stable than S-S-S clusters. S-I-S clusters are found to be strongly bound with respect to their constituents and exhibit partial or total hole compensation. While isolated Mni behaves like a hole-killer and is expected to destroy ferromagnetism, in (MnGa-Mni-MnGa)0, the Mni is found to mediate the ferromagnetic arrangement of spins on MnGa. The charged complex (MnGa-Mni-MnGa)2+ has a similar ferromagnetic stabilization energy on the two MnGa sites as in Mn$`{}_{}{}^{0}{}_{Ga}{}^{}`$-Mn$`{}_{}{}^{0}{}_{Ga}{}^{}`$ cluster without Mni almost as if Mni did not exist. Thus ferromagnetism in Mn doped GaAs arises from holes due to substitutional MnGa, as well as from MnGa-Mni-MnGa complexes.

This work was supported by the U.S. DOE, Office of Science, BES-DMS under contract no. DE-AC36-99-G010337.


format: text

text:
# High pressure phases in highly piezoelectric PbZr0.52Ti0.48O3

\[

## Abstract

Two novel room-temperature phase transitions are observed, via synchrotron x-ray diffraction and Raman spectroscopy, in the PbZr0.52Ti0.48O3 alloy under hydrostatic pressures up to 16 GPa. A monoclinic (M)-to-rhombohedral (R1) phase transition takes place around 2-3 GPa, while this R1 phase transforms into another rhombohedral phase, R2, at $``$ 6-7 GPa. First-principles calculations assign the R3m and R3c symmetry to R1 and R2, respectively, and reveal that R2 acts as a pressure-induced structural bridge between the polar R3m and a predicted antiferrodistortive R$`\overline{3}`$c phase.

\]

PbZr1-xTixO3 (also called PZT) solid solutions with x $``$ 0.5 are used, already for decades, in virtually all piezoelectric devices from ultrasound generators to micropositioners, due to their outstanding electromechanical performance . All PZT compositions have the cubic perovskite structure at high temperatures, while they undergo a phase transition, at about 650 K, to a ferroelectric rhombohedral R3m phase (with a polarization pointing along the  pseudo-cubic direction) for Zr-rich samples and to a ferroelectric tetragonal P4mm phase (with a polarization lying along ) for Ti-rich samples . Recently, higher resolution synchrotron x-ray experiments at low temperatures in high quality samples have made possible the resolution of the region around the R3m-P4mm boundary, and have revealed the existence of a third ferroelectric phase with a lower monoclinic (M) Cm symmetry within a narrow compositional range . This M phase can be considered as a structural bridge between the R3m and P4mm phases since its polarization continuously rotates in the ($`\overline{1}`$10) plane from the  to the  pseudo-cubic directions, as the Ti content increases . Raman measurements in PZT have also found mode splitting which are consistent with monoclinic symmetry , and first-principles calculations have shown that the polarization rotation is responsible for the large piezoelectric coefficients observed in PZT .

Furthermore, antiphase oxygen rotations, superimposed on polar displacements, have been recently observed within the monoclinic phase of PZT at low temperature . This is consistent with first-principles calculations predicting that oxygen rotations and polar displacements are close in energy . Interestingly, the energetical order between different phases in perovskite materials can be modified by applying an external pressure. This is clearly evidenced by earlier works observing that pressure lowers the Curie temperature of phase transitions induced by softening of zone-center phonons, like the $`\mathrm{\Gamma }_{15}`$ polar modes in PbTiO3  — and thus leads to a cubic paraelectric phase at high enough pressure, as recently reported for PbZr0.52Ti0.48O3  — while it increases the temperature of the transitions related to zone-boundary modes like the antiferroelectric mode in PbZrO3 . These opposite effects , the fact that Cm is a structure bridging the R3m and P4mm phases (and is thus very sensitive to a small change of interactions ), and the observed coexistence of polar displacements and oxygen rotation suggest that the phase diagram of PZT versus pressure can be extraordinary rich.

Motivated to precisely determine and fully understand this phase diagram, we combined three complementary techniques, synchrotron x-ray diffraction, Raman spectroscopy and first principles calculations, to investigate PbZr0.52Ti0.48O3 under hydrostatic pressure. These studies revealed (1) two novel phase transitions, (2) two different pressure-induced bridging structures, and (3) the absence of a high pressure cubic paraelectric phase.

The same PbZr0.52Ti0.48O3 ceramic pellets of Ref.  were used. A piece was crushed into powder and loaded in a diamond anvil cell using both nitrogen and a 4:1 methanol-ethanol-mixture as pressure media for the x-ray and Raman measurements, respectively. A maximum pressure of 16 GPa was reached in hydrostatic conditions. The pressure was measured by the ruby fluorescence method . X-ray diffraction data were taken at room temperature in the ID9 beamline at the ESRF with a monochromatic beam of $`\lambda `$ =0.414 Å focused to a 30x30 $`\mu m^2`$ spot and using an angle dispersive set-up with image plate detector. The unit cells for the different pressures were obtained by means of a full profile Le Bail analysis , using a pseudo-voigt peak shape with asymmetry correction . A second cubic phase was used to account for the diffuse scattering found in the monoclinic phase, as in ref. . Depolarized Raman spectra of the sample powders were recorded at room temperature in back-scattering geometry with a Jobin Yvon T64000 spectrometer equipped with a microscope objective and the 514.5 nm line of an Ar+ ion laser as excitation line. The Raman spectra after pressure release were identical to the initial spectra attesting the reversibility of pressure-induced changes.

As shown in Figure 1(top), at atmospheric pressure PbZr0.520Ti0.48O3 is just at the monoclinic-tetragonal phase transition showing a very small monoclinic distortion, with lattice constants am, bm, cm and a $`\beta `$ angle. Pressure induces fairly large and continuous changes in $`c_m`$, in agreement with ref., and far more subtle changes in am, bm and $`\beta `$. Small pressures below 2 GPa increase the monoclinic distortion by increasing both $`\beta `$ and the difference between am and bm. At $``$ 2-3 GPa, am= bm=cm, and the triplet of the pseudo-cubic (222) reflection (see Fig. 1 (bottom) ) observed at low pressures changes into a doublet. The corresponding phase, denoted by R1, is rhombohedral with an angle $`\alpha `$ $`<`$ 90o , or unit cell elongated along  . The most plausible space group for R1 is R3m since this latter is ferroelectric with polarization along . For P$`>`$ 3 GPa, the rhombohedral distortion of R1 decreases with increasing pressure, until at about 5 GPa when the diffraction pattern seems to be cubic with no visible peak split in agreement with ref.. However, an analysis of the peak widths shows that they are not resolution limited, suggesting local non-cubic distortions. Interestingly, when the pressure is further increased, above 7 GPa, the distortion from cubic continuously increases and the diffraction peaks start splitting again. As seen in Figure 1(bottom), the diffraction pattern is again that of a rhombohedral phase but this time with $`\alpha `$$`>`$ 90o . This new phase (denoted by R2) is thus associated to a unit cell flattened along . One good candidate for the symmetry of R2 is R$`\overline{3}`$c since all known rhombohedral perovskites with $`\alpha `$$`>`$ 90o adopt this space group, a well-known example being LaAlO3. In this non-polar structure the distortion from cubic is due to anti-phase rotations of the oxygen octahedra around the  direction, which results from the condensation of zone-boundary modes (R25) and thus double the unit cell . This type of rotations is also observed in Zr-rich PZT at atmospheric pressure and low temperature, although in this case polar displacements are also present and the structure has the space group R3c with $`\alpha `$$`<`$90o .

Raman measurements were performed to complement the x-ray data. Figure 2 presents the evolution of the Raman spectra of the same sample with increasing pressure until 10.2 GPa. A particularly striking finding is that the Raman signal is fairly well-defined for any pressure, including the ones around 5-7 GPa for which x-ray shows a cubic-like structure, implying the existence of non-cubic distortions on the Raman characteristic length scale, too small for x-rays to resolve it. Furthermore, the spectra at low pressures are similar to those reported in the recent literature  for similar compositions. Several changes are noticed upon increasing pressure. The most obvious one being the appearance of several new bands at high pressures (marked with arrows in the figure) and the suppression of some others. Most of the changes appear in between 5.5 and 7.2 GPa, in agreement with the critical pressure observed by x-ray diffraction for the transition to the R2 phase. The strong and sharp band at 380 cm-1 has been observed under pressure in other perovskites and has been associated with rhombohedral symmetry . Its narrow width at high pressures is indicative of the long-range character of the distortion, and the increase in the number of bands is also consistent with a doubling of the unit cell. Furthermore, we observe a softening (inset Figure 2) and progressive disappearing of the two modes in the 200-to-350 cm-1 range, clearly supporting the occurence of a phase transition, most possibly related to fundamental changes of the Zr/Ti-cation displacement.

In order to gain further insights, we also perform 0K calculations within the local density-approximation (LDA) , using the virtual crystal approximation (VCA) proposed in Ref.  (in which further details of the calculations can be found). We limit ourselves to the study of 7 phases: the paraelectric cubic Pm-3m state; the ferroelectric monoclinic Cm, rhombohedral R3m and tetragonal P4mm phases; the antiferrodistortive rhombohedral R$`\overline{3}`$c and tetragonal I4/mcm phases; and the rhombohedral R3c state that can exhibit both polar and antiferrodistortive degrees of freedom. The lattice vectors of the Cm, R3m and P4mm phases are those experimentally determined in Refs.  for an atmospheric pressure, and the associated strain tensors – with respect to the cubic case – are kept frozen in our calculations. On the other hand, for computational convenience, the lattice vectors of the R$`\overline{3}`$c, I4/mcm and R3c phases do not contain any strain with respect to the cubic phase. Fig 3a shows the $`\mathrm{\Delta }H=HH_{cubic}`$ difference between the enthalpy of each phase and the $`H_{cubic}`$ enthalpy of the cubic paraelectric phase as a function of pressure. It is well known that LDA usually underestimates the lattice constant, which explains why the LDA-predicted pressure corresponding to the experimental volume determined in Ref  for the Cm phase, at low temperature and at an atmospherique pressure, is found to be $``$8 GPa. Therefore, in order to compare with the experiments, we have performed an uniform shift of the theoretical $`P`$ scale by 8 GPa, as indicated in Figs 3 by the rescaled theoretical pressure $`\overline{P}`$.

One can notice that our simulations yield a R3m phase that is slightly more favorable than Cm for $`\overline{P}=0`$. This discrepancy between theory and experiments for the ground-state symmetry is likely due to the fact that we neglect an important feature of monoclinic phases in our VCA calculations, namely the significant disorder existing between local dipoles located in different 5-atom cells . On the other hand, further using this rescaled pressure scale leads to a R3m–R3c transition occurring around 6 GPa. This is in remarkable agreement with the experimentally-determined R1–R2 transition pressures mentioned above , and thus strongly suggests that the spaces groups of R1 and R2 are R3m and R3c, respectively. The calculations further predict that R3c transforms into a R$`\overline{3}`$c phase for pressure $`\overline{P}`$ around 20 GPa. (Note that $`\mathrm{\Delta }H`$ of the R3c and R$`\overline{3}`$c phases are quite negative even for $`\overline{P}`$ larger than 20 GPa, which further confirms the absence of a high pressure cubic phase.)

Fig 3b and Fig 3c show the pressure behavior of the local soft-mode $`𝐮`$ (directly related to the spontaneous polarization ) and the rotation angle of the oxygen octahedra, respectively, for the relevant phases. One can clearly see that the z-coordinate of $`𝐮`$ in the Cm phase significantly decreases when increasing $`\overline{P}`$ from 0 to 2 GPa, while the corresponding x- and y-components are barely affected by this change of pressure. (Note that the x-, y- and z- axes are chosen along the pseudo-cubic ,  and  directions, respectively). These simulated behaviors are consistent with the pressure changes of am, bm and cm depicted in Fig 1 for the M phase, since ferroelectrics are well-known to exhibit a coupling between polarization (related to $`𝐮`$) and strain (related to lattice vectors) . Similarly, the continuous decrease predicted for every component of $`𝐮`$ in the R3m phase – when increasing $`\overline{P}`$ from 2 to 6 GPa – is consistent with the experimentally observed increase (see, Fig 1) towards $`90^o`$ for the rhombohedral angle of R1 in this range. Figs 3b and 3c further indicate that in the R3c phase for rescaled pressures above 3 GPa, the spontaneous polarization decreases in magnitude while the rotation of the oxygen octahedra increases, at increasing pressure. First-principles calculations thus provide a successful explanation for the unusual phase that is Raman-active but characterized by a $`\alpha `$ angle $`90^o`$ between its lattice vectors, and that is observed for pressure ranging between 5 and 7 GPa (see Figs 1 and 2): this phase is simply a R3c structure in which polar displacements (favoring $`\alpha <90^o`$ ) compete with the rotation of the oxygen octahedra (favoring $`\alpha >90^o`$). When further increasing pressure, the simulations yield a significant increase (decrease) of the antiferrodistortive (polar) structural features. This explains why $`\alpha `$ increases with pressure after 7 GPa (see Fig.1), and also indicates that R3c behaves as a pressure-induced bridging structure between the R3m and R$`\overline{3}`$c phases (which are not in a simple super- or subgroup relation).

In summary, complementary synchrotron x-ray diffraction, Raman spectroscopy and first principles calculations in piezoelectric PbZr0.52Ti0.48O3, under pressure, reveal two novel room-temperature phase transitions. The first transition occurs around 2-3 GPa from the ferroelectric monoclinic Cm phase to a ferroelectric rhombohedral R3m phase. The second transition occurs at $``$6-7 GPa, from the ferroelectric rhombohedral R3m phase to a rhombohedral R3c phase exhibiting both polar and antiferrodistortive structural features. The first transition thus involves polarization rotation, while the second transition leads to a smooth suppression of the polar displacements and continuously increases the oxygen rotations around the Zr/Ti cations with increasing pressure. The pressure phase diagram of PZT thus contains two peculiar bridging structures: polar Cm that can be seen as the natural evolution of polar P4mm to polar R3m, and R3c that brings polar R3m to antiferrodistortive R$`\overline{3}`$c. The absence of a high-pressure cubic phase is, therefore, reported. Our results may also be relevant to thin films because of the crucial role played by the stress in these two-dimensional systems. Finally, we would like to raise the possibility that the mysterious (cubic-like while polar) so-called X-phase, recently reported in some complex perovskites , is a phase similar to the one we observe at the R3m-R3c transition, namely a phase in which polar distortions compete with antiphase oxygen rotations.

We would like to thank R. Guo and the late S.K. Eagle Park for the excellent samples and D.E. Cox and G. Shirane for very useful discussions. The theoretical part of this work is supported by Office of Naval Research Grants N00014-01-1-0600 and N00014-01-1-0365 (Center for Piezoelectrics by Design), and National Science Foundation Grant DMR-9983678.