text:
# First-principles calculation of the effect of strain on the diffusion of Ge adatoms on Si and Ge (001) surfaces

\]http://www.mit.edu/ avdw/

## Abstract

First-principles calculations are used to calculate the strain dependencies of the binding and diffusion-activation energies for Ge adatoms on both Si(001) and Ge(001) surfaces. Our calculations reveal that the binding and activation energies on a strained Ge(001) surface increase and decrease, respectively, by 0.21 eV and 0.12 eV per percent compressive strain. For a growth temperature of 600 C, these strain-dependencies give rise to a 16-fold increase in adatom density and a 5-fold decrease in adatom diffusivity in the region of compressive strain surrounding a Ge island with a characteristic size of 10 nm.

In heteroepitaxial growth of lattice-mismatched thin films, the Stranski-Krastanov (SK) growth mode has been widely investigated as a basis for self-assembling arrays of coherent nanostructured islands, commonly referred to as quantum dots (QDs).Eaglesham and Cerullo (1990) The need for highly monodisperse QD arrays in semiconductor optoelectronic device applications has motivated extensive research into the microscopic mechanisms influencing the evolution of island size distributions during SK growth. While the stabilizing thermodynamic effects associated with the elastic interactions between strained islands have been investigated in some detail (e.g., Refs. Shchukin and Bimbert, 1999; Daruka and Barabási, 1997, and references cited therein), the role of various proposed kinetic mechanisms on the development of island size distributions remains less clear.

Within a self-consistent mean-field rate theory, Koduvely and Zangwill Koduvely and Zangwill (1999) demonstrated that with decreasing island-island separation, a strain-mediated decrease in the barrier for adatom-island detachment leads to a reduction in the mean island size, with an associated narrowing of the size distribution. These findings are qualitatively consistent with experimental observations Kobayashi et al. (1996) in InAs/GaAs. Madhukar Madhukar (1996) and Penev, Kratzer and Scheffler Kratzer et al. (2002); Penev et al. (2001) have considered the growth of InAs QDs on GaAs, where increasing island size leads to a buildup of compressive elastic strains in the surrounding substrate.Moll et al. (1998) Within simplified models of diffusion-limited growth these authors demonstrated that, in the InAs/GaAs heteroepitaxial system, the strain-dependence of the parameters governing adatom diffusion gives rise to a reduction in the flux of adatoms reaching larger islands relative to small ones, leading to a reduced rate of coarsening and an associated narrowing of the size distribution.

While the potentially important consequences for island growth kinetics arising from strain dependencies in adatom binding and migration energies have been clearly demonstrated, attempts to determine the magnitude of these effects in specific systems have been undertaken in relatively few systems.Penev et al. (2001); Ratsch et al. (1997); Shu et al. (2001); Roland and Gilmer (1992); Spjut and Faux (1994); Schroeder and Wolf (1997) These effects are most readily investigated using first-principles calculations Kaxiras (1996); Smith et al. (1995); Penev et al. (2001); Ratsch et al. (1997); Shu et al. (2001), as they are extremely difficult to isolate experimentally. The purpose of the present work is to investigate the effect of strain-dependent adatom diffusion and binding energies upon island growth kinetics in Ge/Si(001). This system represents one of the most widely studied examples of QD formation induced by SK growth, Ross et al. (1998, 1999); Williams et al. (2000); Floro et al. (1997); Mo et al. (1990); Tersoff et al. (1996) yet to date the effect of strain upon Ge adatom binding energies and diffusion rates on Ge wetting layers remains unstudied. We employ first-principles calculations to compute the strain dependence of Ge adatom energetics on Si(001) and Ge(001) surfaces.

Our calculations were set up as follows. We consider adatom diffusion in the dilute limit (as opposed to dimer diffusion Zoethout et al. (2000)). The initial atomic configuration for all of our calculations is taken to be the well-known minimum energy $`c(4\times 2)`$ reconstruction (depicted in the upper left corner of Fig. 1). The adatom calculations are performed using a supercell geometry where an artificial three-dimensional periodicity is imposed on the system. This artificial periodicity is harmless, provided that convergence with respect to the distance between the periodic images is achieved. Following earlier studies of adatom diffusion of similar systems, Milman et al. (1994, 1996); Milman (1997) our supercell consists of two repetitions of the surface unit cell, depicted in Fig. 1 along the dimer rows direction (i.e. a $`4\times 4`$ supercell). Twelve layers of atoms were used, separated by 6 layers of vacuum. By varying the thickness of the slab and the width of the vacuum separating the periodic images perpendicular to the surface, it was verified that the selected system size provides energies with an accuracy of the order of 0.03 eV. Single adatoms were placed on each of the two slab surfaces, which has the effect of doubling the adatom’s energetic contribution, thereby improving the accuracy of the calculations. All calculations were performed using the ab initio total-energy and molecular-dynamics program VASP (Vienna ab initio simulation package) developed at the Institut für Materialphysik of the Universität Wien Kresse and Furthmüller (1996a, b), which implements Vanderbilt ultra-soft Vanderbilt (1990) pseudo-potentials Phillips and Kleinman (1959). The energy cutoff for the plane-wave basis set was set to 150 eV and a $`2\times 2\times 1`$ mesh of $`k`$-points was used for the $`4\times 4`$ supercell. All atoms were allowed to relax.

On the Si $`(001)`$ surface, the diffusion of a Ge adatom along the dimer rows is known to be of the order of 1000 times faster than across dimer rows Mo and Lagally (1991) and we therefore focus solely on the diffusion along dimer rows. Since a typical surface consists of terraces where the direction of the dimer rows changes by 90 at each monoatomic step, fast diffusion in any direction is possible at the mesoscopic level, even though the fast diffusion is unidirectional at the microscopic level. Earlier studies Milman et al. (1994, 1996); Milman (1997); Dalpian et al. (2001) have unambiguously determined the location of the binding site and the activated state of the Ge adatom on the Si $`(001)`$ surface (see Fig. 1) as well as the precise configuration of the surface dimers in the vicinity of the adatom.

Since experimental and computational evidence is scarcer in the case of the Ge adatom on the Ge $`(001)`$ surface, Mea (2001) we considered 7 possible binding sites and found the minimum energy site to be as shown in Figures 1 and 2. This site remains energetically favored for values of the substrate lattice parameter up to 3% larger than the Si lattice parameter. Beyond that threshold, a binding site analogous to the binding site of Ge on Si (001) becomes favorable. We shall neglect this alternate binding site, as the substrate lattice parameter required to stabilize it falls outside the range sampled during Ge on Si (001) heteroepitaxy.

We identified the diffusion path using the nudged elastic band method G. Mills and Schenter (1995); H. Jonsson and Jacobsen (1998) and refined the position of the saddle point (see Figures 1 and 2) using the quasi-Newton algorithm. Kresse and Furthmüller (1996b) We focused on the diffusion along the valleys between the dimer rows, because the binding sites located atop the dimer row had a binding energy at least 0.3 eV larger than the saddle point energy for diffusion within the valleys, strongly suggesting that any other saddle points located on the dimer rows would also have a higher energy than the saddle point we analyzed. The nearly one-dimensional diffusion and the adatom’s preference for sites located in the valley between the dimer rows agrees qualitatively with the results of earlier calculations based on semi-empirical potentials, Mea (2001) although the precise location of the binding sites and of the saddle points differ.

For each type of surface (Si or Ge), the energies of three geometries were calculated (the free surface, the saddle point and the binding site configurations) at various levels of biaxial strain imposed parallel to the plane of the surface. The resulting energy versus strain relationships are plotted in Fig. 3. Our results reveal four important observations. First, linear approximations to the strain-dependence of the binding energies and activation barriers Shu et al. (2001) must clearly be used with care. Second, the Ge adatom binding energy does not necessarily exhibit a minimum when the lattice parameter of the substrate matches the lattice parameter of bulk Ge. Third, the sign of the change of these energies under strain are highly system-dependent: calculations on the structurally similar In/GaAs(001) heteroepitaxial system Penev et al. (2001) found the strain-dependences of the binding ($`E_b`$) and saddle-point energies ($`E_b+E_a`$) the be of a sign opposite to the ones in the Ge/Ge(001) system. The sign of the strain-dependence of the activation barrier ($`E_a`$) for Ge on both Si (001) and Ge (001), however, agrees with earlier studies in In/GaAs(001) Penev et al. (2001) and Si/Si(001) Shu et al. (2001); Roland and Gilmer (1992); Spjut and Faux (1994); Zoethout et al. (2000). Finally, the magnitude of the effect of strain on binding and saddle point energies found in the case of Ge on Ge(001) is by far the largest, relative to any other system studied so far.Penev et al. (2001); Shu et al. (2001); Ratsch et al. (1997); Roland and Gilmer (1992); Spjut and Faux (1994); Zoethout et al. (2000)

In order to quantify the importance of these results in the context of quantum dot growth, we calculated the strain field in the vicinity of a Ge $`(105)`$-terminated pyramidal island with a characteristic size of 10 nm.Medeiros-Ribeiro et al. (1998) For the calculation of strain fields around coherent strained islands, both finite-element techniques (e.g., Ref. Moll et al., 1998) as well as approximate analytical methods Tersoff and Tromp (1993) have been applied previously. In the present work, elastic strain fields have been derived from the results of the linear stability analysis of Spencer, Voorhees and Davis.Spencer et al. (1993) In this analysis, the linearized strain fields arising from Fourier-mode perturbations in the surface height were derived within isotropic elasticity theory. By summing the resulting expressions over the Fourier amplitudes describing the shape function of a Ge pyramid, we calculated the in-plane strain field (shown in Fig. 4) at the surface of a Ge wetting layer (three atomic-planes in thickness) on Si(001). In Fig. 4 the strain is referenced to the Si lattice constant, so that zero corresponds to a wetting layer epitaxially strained on the Si substrate. In agreement with previous calculations for related systems, Moll et al. (1998); Penev et al. (2001); Tersoff and Tromp (1993) the strain field surrounding the island is compressive in nature.

In Fig. 4, the strain is seen to range in magnitude from zero to roughly $`1`$% in the vicinity of the Ge(105)-terminated island. Over this range of compressive strains, the binding energy of a Ge adatom on a Ge surface varies by about 0.21 eV, while the activation barrier varies by about 0.12 eV, as seen in Fig. 3. At a typical deposition temperature of 600 C, these variations correspond, respectively, to a 16-fold variation in equilibrium adatom density and a 5-fold variation in adatom mobility. In contrast, the corresponding strain-dependences on the Si (001) surface are much less pronounced. Of course, the strain-dependence of the entropic prefactors Penev et al. (2001) could very well modify those order-of-magnitude estimates based solely on the strain-dependence of the energetic contributions.

To further explore the implications of our calculated results, we follow the analysis of Penev, Kratzer and Scheffler, Penev et al. (2001) employing a simple surface diffusion-limited one-dimensional model of adatom diffusion between two neighboring islands. We expand slightly on their work by including boundary conditions at the island edges that impose local equilibrium between adatoms and islands. For simplicity, we consider all entropic prefactors to be strain-independent. We perform a stability analysis by considering two islands of equal size located at $`x=0`$ and $`x=l`$ and calculate the changes in the adatom flux toward each island as the size of one is perturbed. Let $`F_0,F_l`$ denote the flux towards each island and consider the (dimensionless) asymmetry in the flux towards each island, $`F=\left(F_lF_0\right)/\left(2\varphi l\right)`$, where $`\varphi `$ is the adatom deposition rate and $`l`$ is the distance separating the two islands. Solving the above diffusion problem leads to

$$F=\frac{K}{M_0\varphi l}\frac{\left(e^{\beta E_{i,l}}e^{\beta E_{i,0}}\right)}{l}\frac{M_1}{M_0},$$
(1)
where

$$M_n=\frac{1}{l^{n+1}}_0^l\left(x\frac{l}{2}\right)^ne^{\beta \left(E_a\left(x\right)+E_b\left(x\right)\right)}𝑑x.$$
(2)
and where $`E_a\left(x\right)`$ and $`E_b\left(x\right)`$ are, respectively, the activation barrier and binding energy as a function of position $`x`$; $`E_{i,0}`$ and $`E_{i,l}`$ are, respectively, the energies of an atom bound to the islands located at $`x=0`$ and $`x=l`$; $`K`$ is a constant incorporating all entropic prefactors and $`\beta `$ is reciprocal temperature $`\left(k_BT\right)^1`$.

Our analysis focuses on the diffusion of Ge adatoms on the Ge (001) surface, strained to match the lattice constant of the Si substrate, since a Ge wetting layer covering the Si substrate is known to form in SK growth. We assume that adatoms are insensitive to the presence of Si under the Ge wetting layer, that no substantial Si-Ge interdiffusion occurs and that surface reconstructions are unaffected by epitaxial strain. Accounting for the latter (e.g., through the inclusion of a strain-dependent concentration of “missing dimers”Liu et al. (1997)) may provide another source of strain-dependent diffusion deserving further consideration.

Under the above assumptions, the quantity $`E_a\left(x\right)+E_b\left(x\right)`$ decreases under a compressive strain, as shown in the middle panels of Fig. 3. When the two islands are identical, $`F=0`$, since $`E_{i,l}=E_{i,0}`$ and $`M_1=0`$ (as $`E_a\left(x\right)+E_b\left(x\right)`$ is then symmetric with respect to $`x=\frac{l}{2}`$). Now, consider an increase in the size of the island at $`x=l`$. The first term of Equation (1) describes the thermodynamic driving force for coarsening: $`e^{\beta E_{i,l}}`$ decreases due to capillarity (a larger island has a smaller surface to volume ratio). This term thus causes $`F`$ to become positive, increasing the flux towards the larger island. This term is inversely proportional to $`\varphi `$, demonstrating that increasing the deposition rate would act to reduce this natural coarsening effect.

The second term of Equation (1), which is independent of deposition flux, quantifies the effect of strain-dependent adatom binding and migration energies. A slight increase in the size of the island at $`x=l`$ induces a compressive strain in the wetting layer in the vicinity of that island, thus inducing a decrease in the saddle point energy $`E_a(x)+E_b(x)`$ in the same area and causing $`M_1`$ to become negative. The result is a relative increase in the adatom flux towards the larger island. In other words, this analysis suggests that the strain-dependencies of Ge adatom binding and migration energies plotted in Fig. 3 act to accelerate the rate of coarsening of larger islands.

In summary, first-principles calculations have been employed to compute the strain dependence of Ge adatom binding energies and activation energies for diffusion on both Si(001) and Ge(001) surfaces. For a growth temperature of 600 C, these strain-dependencies give rise to a 16-fold increase in adatom density and a 5-fold decrease in adatom diffusivity in the region of compressive strain surrounding a Ge island with a characteristic size of 10 nm. Within a simplified model of diffusion-limited growth, these strain dependencies are found to have the qualitative effect of accelerating the natural coarsening rate of larger islands, in contrast with earlier results obtained for the related InAs on GaAs system.Penev et al. (2001) The large magnitude of the strain dependence of adatom energetics on Ge(001) obtained in the present calculations also contrasts with earlier findings in other systemsRatsch et al. (1997); Roland and Gilmer (1992); Spjut and Faux (1994); Penev et al. (2001); Shu et al. (2001); Zoethout et al. (2000) and warrants further consideration of these effects in more detailed kinetic models to further elucidate their effect on island growth.

###### Acknowledgements.

This work was supported by the NSF under programs DMR-0102794 and NSF-MRSEC DMR-00706097, using computer resources provided by the National Partnership for Advanced Computational Infrastructure at the University of Michigan.


format: text

text:
# References

Anomalous net-baryon-rapidity spectra at RHIC
Georg Wolschin111E-mail: wolschin@uni-hd.de http://wolschin.uni-hd.de
Institut für theoretische Physik der Universität, D-69120 Heidelberg, Germany

(Submitted January 02, 2003)

Net-baryon rapidity distributions in central Au+Au collisions at $`\sqrt{s_{NN}}`$ = 200 GeV are shown to consist of three components. The nonequilibrium contributions are accounted for in a Relativistic Diffusion Model. Near midrapidity, a third fraction containing $`Z_{eq}22`$ protons reaches local statistical equilibrium in a discontinuous transition. It may be associated with a deconfinement of the participant partons and thus, serve as a signature for Quark-Gluon Matter formation.

Rapidity distributions of participant (net) baryons are very sensitive to the dynamical and statistical properties of nucleus-nucleus collisions at high energies. Recent results for net-proton rapidity spectra in central Au+Au collisions at the highest RHIC energy of $`\sqrt{s_{NN}}`$ = 200 GeV show an unexpectedly large rapidity density at midrapidity. The BRAHMS collaboration finds  dN/dy = 7.1 $`\pm `$ 0.7 (stat.) $`\pm `$ 1.1 (sys.) at y = 0.

The $`\mathrm{\Lambda },\overline{\mathrm{\Lambda }}`$ feed-down corrections reduce this yield by 17.5 per cent  when performed in accordance with the PHENIX $`\mathrm{\Lambda }`$results  at 130 GeV, but the amount of stopping remains significant, although a factor of about 4 smaller as compared to Pb+Pb at the highest SPS energy. (A corresponding STAR result  for y=0 at 130 GeV does not yet include the feed-down correction). Many of the available numerical microscopic models encounter difficulties to predict the net-proton yield in the central midrapidity valley of the distribution, together with the broad peaks at the detected positions.

In this Letter I propose to interpret the data in a nonequilibrium-statistical Relativistic Diffusion Model. The net baryon rapidity distribution at RHIC energies emerges from a superposition of the beam-like nonequilibrium components that are broadened in rapidity space through diffusion due to soft (hadronic, low $`p_{}`$) collisions and particle creations, and a statistical equilibrium (thermal) component at midrapidity that arises from hard (partonic, high $`p_{}`$) processes.

At RHIC energies, the underlying distribution functions turn out to be fairly well separated in rapidity space. Since the transverse degrees of freedom are in (or very close to) thermal equilibrium, they are expected to decouple from the longitudinal ones. The time evolution of the distribution functions is then governed by a Fokker-Planck equation (FPE) in rapidity space . In the more general case of nonextensive (non-additive) statistics  that accounts for long-range interactions and violations of Boltzmann’s Stosszahlansatz  as well as for non-Markovian memory (strong coupling) effects , the resulting FPE for the rapidity y in the center-of-mass frame is

$$\frac{}{t}[R(y,t)]^\mu =\frac{}{y}\left[J(y)[R(y,t)]^\mu \right]+D(t)\frac{^2}{y^2}[R(y,t)]^\nu .$$
(1)
Since the norm of the rapidity distribution has to be conserved, $`\mu `$ = 1 is required, and the nonextensivity parameter that governs the shape of the power-law equilibrium distribution becomes q = 2 - $`\nu `$ . In statistical equilibrium, transverse mass spectra and transverse momentum fluctuations in relativistic systems at SPS-energies $`\sqrt{s_{NN}}`$=17.3 GeV require values of q very slightly above one, typically q = 1.038 for produced pions in Pb+Pb . For $`q1`$, the equilibrium distribution converges to the exponential Boltzmann form, whereas for larger values of q (with $`q<1.5`$) significantly broader equilibrium distributions are obtained, and the time evolution towards them becomes superdiffusive .

To study rapidity distributions in multiparticle systems at RHIC energies in a nonequilibrium-statistical framework , I start with q = $`\nu `$ = 1 corresponding to the standard FPE. For a linear drift function

$$J(y)=(y_{eq}y)/\tau _y$$
(2)
with the rapidity relaxation time $`\tau _y`$, this is the so-called Uhlenbeck-Ornstein process, applied to the relativistic invariant rapidity. The equilibrium value is $`y_{eq}=0`$ in the center-of-mass for symmetric systems, whereas $`y_{eq}`$ is calculated from the given masses and momenta for asymmetric systems. Using $`\delta `$function initial conditions at the beam rapidities $`\pm y_b`$ ($`\pm `$5.36 at p=100 GeV/c per nucleon), the equation has analytical Gaussian solutions. The mean values shift in time towards the equilibrium value according to

$$<y_{1,2}(t)>=y_{eq}[1exp(2t/\tau _y)]\pm y_b\mathrm{exp}(t/\tau _y).$$
(3)
For a constant diffusion coefficient $`D_y`$, the variances of both distributions have the well-known simple form

$$\sigma _{1,2}^2(t)=D_y\tau _y[1\mathrm{exp}(2t/\tau _y)],$$
(4)
whereas for a time dependent diffusion coefficient $`D_y(t)`$ that accounts for collective (multiparticle) and memory effects the analytical expression for the variances becomes more involved . At short times $`t/\tau _y<<1`$, a statistical description is of limited validity due to the small number of interactions. A kinematical cutoff prevents the diffusion into the unphysical region $`|y|>y_b`$. For larger values of $`t/\tau _y`$, the system comes closer to statistical equilibrium such that the FPE is valid.

Since the equation is linear, a superposition of the distribution functions emerging from $`R_{1,2}(y,t=0)=\delta (yy_b)`$ yields the exact solution, with the normalization given by the total number of net baryons and the value of $`t/\tau _y`$ at the interaction time t=$`\tau _{int}`$ (the final time in the integration of (1)) determined by the peak positions . This approach has also been applied sucessfully to produced particles at RHIC energies , although there the initial conditions are less straightforward.

The microscopic physics is contained in the diffusion coefficient. Macroscopically, the transport coefficients are related to each other through the dissipation-fluctuation theorem (Einstein relation) with the equilibrium temperature T

$$D_y=\alpha Tf(\tau _y,T).$$
(5)
In  I have obtained the analytical result for $`D_y`$ as function of T and $`\tau _y`$ from the condition that the stationary solution of (1) is equated with a Gaussian approximation to the thermal equilibrium distribution in y-space (which is not exactly Gaussian, but very close to it) as

$$D_y(\tau _y,T)=\frac{1}{2\pi \tau _y}\left[c(\sqrt{s},T)m^2T(1+2\frac{T}{m}+2(\frac{T}{m})^2)\right]^2exp(\frac{2m}{T})$$
(6)
with $`c(\sqrt{s},T)`$ given in  in closed form. This allows to maintain the linearity of the model and hence, to solve the FPE analytically, although small corrections are to be expected. They cause minor deviations in the calculated rapidity distributions that are within the size of the error bars of the experimental data at SPS energies.

In the linear model, net baryon rapidity spectra at low SIS-energies are well reproduced, whereas at AGS-, SPS- and RHIC energies I find discrepancies to the data that rise strongly with $`\sqrt{s}`$. At SPS energies, this has recently been confirmed in a numerical calculation  based on a nonlinear drift

$$J(y)=\alpha m_{}sinh(y)\alpha p_{}$$
(7)
with the transverse mass $`m_{}=\sqrt{m^2+p_{}^2}`$, and the longitudinal momentum $`p_{}`$. Together with the dissipation-fluctuation theorem (5), this yields exactly the Boltzmann distribution as the stationary solution of (1) for $`\nu =q=1`$. The corresponding numerical solution with $`\delta `$function initial conditions at the beam rapidities is, however, only approximately correct since the superposition principle is not strictly valid for a nonlinear drift. Still, the numerical result shows almost the same large discrepancy between data and theoretical rapidity distribution as the linear model. In a q=1 framework, the net proton distribution in Pb+Pb at the highest SPS energy requires a rapidity width coefficient $`\sqrt{D_y\tau _y}`$ that is enhanced beyond the theoretical value (5) by a factor of $`g(\sqrt{s})2.6`$ due to memory and collective effects , Fig.1.

Alternatively, a transition to nonextensive statistics  maintaining the weak-coupling diffusion coefficient from (5) requires a value of q that is significantly larger than one. In an approximate numerical solution of (1) with the nonlinear drift (7), q=1.25 has been determined for the net-proton rapidity distribution in Pb+Pb collisions at the SPS . The only free parameter is q, whereas in the linear q=1 case the enhancement of $`D_y`$ beyond (5) is the only parameter.

This value of q in the nonlinear model is considerably larger than the result $`q(\sqrt{s_{NN}})`$=1.12 extrapolated from Wilk et al.  at the SPS-energy $`\sqrt{s_{NN}}`$ = 17.3 GeV. Here, the relativistic diffusion approach is applied to produced particles in proton-antiproton collisions in the energy range $`\sqrt{s}`$ = 53 GeV - 1800 GeV, and used to predict LHC-results. The nonlinearity $`q>1`$ appears to be an essential feature of the $`p\overline{p}`$ data. The larger value of q in heavy systems as compared to $`p\overline{p}`$ at the same NN-center-of-mass energy emphasizes the increasing superdiffusive effect of multiparticle collisions both between participants, and between participants and produced particles. It is, however, conceivable that both a violation of (5) due to memory effects, and $`q>1`$ have to be considered in a complete description.

Investigating the Au+Au system at RHIC energies in the linear q=1 model, it is now easy to show that in central collisions (10 per cent of the cross-section) the nonequilibrium net-proton rapidity spectrum calculated with the standard dissipation-fluctuation theorem underpredicts the widths of the nonequilibrium fractions of the experimental distribution  significantly. In this comparison, the temperature T=170 MeV is taken from a thermal fit of charged antiparticle-to-particle ratios in the Au+Au system at 200 GeV per nucleon , and the theoretical value of the rapidity width coefficient calculated from the analytical expression (6) is $`\sqrt{D_y\tau _y}=7.610^2.`$ As was shown in  for SPS results, the discrepancy persists in case of the nonlinear drift (7) that yields the exact Boltzmann-Gibbs equilibrium solution for q=1

$$E\frac{d^3N}{d^3p}=\frac{d^3N}{dym_{}dm_{}d\varphi }Eexp(E/T).$$
(8)
An enhancement of the weak-coupling rapidity width coefficient by a factor of $`g(\sqrt{s})3.7`$ due to collective and memory effects in the system corresponding to a violation of (5) yields a good reproduction of the nonequilibrium contributions with $`\tau _{int}/\tau _y`$=0.26, but the midrapidity valley that is present in the data is completely absent in the extensive nonequilibrium q=1 case, solid curves in Fig.1 (bottom). This remains true in the nonextensive case ($`1<q<1.5`$), with an approximate distribution function  that is given by a linear superposition of Tsallis-like power-law solutions of (1)

$$R_{1,2}(y,t)=[1(1q)\frac{m_{}}{T}cosh(y<y_{1,2}(t)>)]^{\frac{1}{1q}}.$$
(9)
The dashed curves in Fig.1 show the result for q=1.4, T=170 MeV and a mean transverse mass $`<m_{}>`$=1.2 GeV. This solution is far from the nonextensive equilibrium distribution which would be reached for $`<y_{1,2}(t\mathrm{})>=y_{eq}`$, and it is significantly below the midrapidity data. The result is even worse for larger values of $`m_{}`$. In contrast, the Pb+Pb data at SPS energies  are well described both in the linear model  (Fig.1, top) and in the nonlinear case (cf.  for results with a time-dependent temperature and an integration over transverse mass).

It turns out, however, that the RHIC data can be interpreted rather precisely in the linear q=1 framework with the conjecture that a fraction of $`Z_{eq}22`$ net protons near midrapidity reaches local statistical equilibrium in the longitudinal degrees of freedom. The variance of the equilibrium distribution $`R_{eq}(y)`$ at midrapidity is broadened as compared to the Boltzmann result (dashed curve in Fig.2) due to collective (multiparticle) effects by the same factor that enhances the theoretical weak-coupling diffusion coefficient derived from (5). This may correspond to a longitudinal expansion velocity of the locally equilibrated subsystem as accounted for in hydrodynamical descriptions. In the nonextensive model, the corresponding equilibrium distribution is broadened (blue-shifted) according to $`q1.4`$.

Microscopically, the baryon transport over 4-5 units of rapidity to the equilibrated midrapidity region is not only due to hard processes acting on single valence (di)quarks that are described by perturbative QCD, since this yields insufficient stopping. Instead, additional processes such as the nonperturbative gluon junction mechanism  are necessary to produce the observed central valley. This may lead to substantial stopping even at LHC energies where the separation of nonequilibrium and equilibrium net baryon fractions in rapidity space is expected to be even better than at RHIC. In the late thermalization stage , nonperturbative approaches to QCD thermodynamics are expected to be important.

Recent work indicates that one may account for the observed stopping in heavy-ion collisions at SPS and RHIC energies with string-model parameters determined from hadron-hadron collisions . If this was confirmed, the corresponding rapidity distributions would not be considered to be anomalous from a microscopic point of view. However, this view does not offer a distinction between nonequilibrium and equilibrium contributions to the net baryon rapidity spectra, which both exist at RHIC energies, and are anomalously broadened.

Macroscopically, the complete solution of (1) in the q=1 case is a linear superposition of nonequilibrium and equilibrium distributions (Fig.2, bottom)

$$R(y,t=\tau _{int})=R_1(y,\tau _{int})+R_2(y,\tau _{int})+R_{eq}(y).$$
(10)
It yields a good representation of the preliminary BRAHMS data . (In the $`q>1`$ case, the corresponding solution is questionable because the superposition principle is violated). Based on (10), the transition from net-proton rapidity spectra with a central plateau in Pb+Pb at the lower SPS energies , via a double-humped distribution at the maximum SPS energy  to the central valley at RHIC  is well understood. It has not yet been possible to identify a locally equilibrated subsystem of net baryons at midrapidity below RHIC-energies, although it cannot be excluded that it exists. At SPS energies, the data  are well described by the nonequilibrium distributions, and it is much more difficult (and probably impossible) to identify a locally equilibrated component because the relevant rapidity region is comparatively small, and an equilibrated contribution cannot be separated from the nonequilibrium components in rapidity space. In $`p\overline{p}`$-collisions at $`\sqrt{s}=`$53-900 GeV, no convincing signatures of a phase transition were found .

Most remarkably, Fig.2 suggests that in central Au+Au collisions at $`\sqrt{s_{NN}}`$= 200 GeV there is no continuous transition from the nonequilibrium to the equilibrium contribution in net-proton rapidity spectra as function of time. This may well be due to a sudden enhancement in the number of degrees of freedom as encountered in the deconfinement of participant partons, which enforces a very rapid local equilibration in a fraction of the system. The central valley in net-proton rapidity spectra at RHIC energies could thus be used as an indicator for partonic processes that lead to a baryon transfer over more than 4 units of rapidity, and for quark-gluon plasma formation.

To conclude, I have interpreted recent results for central Au+Au collisions at RHIC energies in a Relativistic Diffusion Model (RDM) for multiparticle interactions based on the interplay of nonequilibrium and equilibrium (”thermal”) solutions. In the linear version of the model, analytical results for the rapidity distribution of net protons in central collisions have been obtained. The anomalous enhancement of the diffusion in rapidity space as compared to the expectation from the weak-coupling dissipation-fluctuation theorem due to high-energy multiparticle effects has been discussed using extensive and nonextensive statistics.

A significant fraction of about 14 per cent of the net protons reaches local statistical equilibrium in a fast and discontinuous transition which is likely to indicate parton deconfinement. The precise amount of protons in equilibrium is related to the experimental value of the rapidity density close to y=0 and hence, possible changes in the final data will affect the percentage. It has not yet been possible to isolate a corresponding fraction of longitudinally equilibrated net protons in the Pb+Pb system at SPS energies. Since no signatures of a transition to the quark-gluon plasma have been observed in $`p\overline{p}`$-collisions, Quark-Matter formation is clearly a genuine many-body effect occuring only in heavy systems at sufficiently high energy density. Consequently, a detailed investigation of the flat midrapidity valley found at RHIC, and of its energy dependence is very promising.

Figure captions

Nonequilibrium contributions to the net-proton rapidity spectra of Au+Au at $`\sqrt{s_{NN}}`$ = 200 GeV in the Relativistic Diffusion Model (RDM) with an equilibrium temperature of T=170 MeV (bottom). The solid curve is obtained in the linear model (q=1) with an enhanced width  due to multiparticle effects according to the preliminary BRAHMS data , squares. The dashed curve corresponds to q=1.4 and $`<m_{}>`$=1.2 GeV in the nonlinear model. At SPS energies (top), NA49 data  for central events (5 per cent) including $`\mathrm{\Lambda }`$ feed-down corrections are compared with the linear model .

Net-proton rapidity spectra for central collisions of Au+Au at $`\sqrt{s_{NN}}`$ = 200 GeV consist of two nonequilibrium components (solid peaks, top) plus an equilibrium contribution at T=170 MeV, dashed curve. It is broadened due to collective (multiparticle) effects, shaded area, and after hadronization, it contains $`Z_{eq}22`$ protons. Superposition creates the flat valley near midrapidity (bottom) in agreement with the preliminary BRAHMS data points ; diamonds include $`\mathrm{\Lambda }`$ feed-down corrections at y=0 (17.5 per cent) and y=2.9 (20 per cent), respectively. Arrows indicate the beam rapidities $`\pm y_b`$.


text:
# Spin-independent origin of the strongly enhanced effective mass in a dilute 2D electron system

## Abstract

We have accurately measured the effective mass in a dilute two-dimensional electron system in silicon by analyzing temperature dependence of the Shubnikov-de Haas oscillations in the low-temperature limit. A sharp increase of the effective mass with decreasing electron density has been observed. Using tilted magnetic fields, we have found that the enhanced effective mass is independent of the degree of spin polarization, which points to a spin-independent origin of the mass enhancement and is in contradiction with existing theories.

The ground state of an ideal, strongly interacting two-dimensional (2D) electron system is predicted to be Wigner crystal chapl . The strength of the interactions is usually characterized by the ratio between the Coulomb energy and the Fermi energy, $`r_s=E_c/E_F`$. Assuming that the effective electron mass is equal to the band mass, the interaction parameter $`r_s`$ in the single-valley case reduces to the Wigner-Seitz radius, $`1/(\pi n_s)^{1/2}a_B`$ and therefore increases as the electron density, $`n_s`$, decreases (here $`a_B`$ is the Bohr radius in semiconductor). In the strongly-interacting limit, no analytical theory has been developed to date. According to numeric simulations tanatar , Wigner crystallization is expected in a very dilute regime, when $`r_s`$ reaches approximately 35. The refined numeric simulations atta  have predicted that prior to the crystallization, in the range of the interaction parameter $`25r_s35`$, the ground state of the system is a strongly correlated ferromagnetic Fermi liquid. At yet higher electron densities, at $`r_s1`$, the (weakly-interacting) electron liquid is expected to be paramagnetic, with the effective mass, $`m`$, and Landé $`g`$ factor renormalized by interactions. Enhancement of $`g`$ and $`m`$ within the Fermi liquid theory is due to spin exchange effects, with renormalization of the $`g`$ factor being dominant compared to that of the effective mass renorm . In contrast, the dominant increase of $`m`$ near the onset of Wigner crystallization follows from an alternative description of the strongly-interacting electron system beyond the Fermi liquid approach, which also predicts the renormalization of $`m`$ to be strongly sensitive to the polarization of spins spivak ; dol .

In dilute silicon metal-oxide-semiconductor-field-effect-transistors (MOSFETs), a strong metallic temperature dependence of the resistance was observed a decade ago krav . Although this anomaly was almost immediately attributed to strong electron-electron interactions, only after a strongly enhanced ratio of the spin and the cyclotron splittings was found at low $`n_s`$ ssc  has it become clear that the system behaves well beyond the weakly interacting Fermi liquid. Later, it was reported that the magnetic field required to produce complete spin polarization, $`B_cn_s/gm`$, tends to vanish at a finite electron density $`8\times 10^{10}`$ cm-2 sh ; vit . These findings point to a sharp increase of the spin susceptibility and possible ferromagnetic instability in dilute silicon MOSFETs. In very dilute GaAs/AlGaAs heterostructures, a similar behavior has been observed in both the 2D hole gao  and electron zhu  system. Recently, experimental results have indicated that in silicon MOSFETs it is the effective mass, rather than the $`g`$ factor, that sharply increases at low electron densities sh1 . The crucial point for understanding these properties of dilute 2D electron liquids is how the effective mass changes with the degree of spin polarization.

In this Letter, we report accurate measurements of the effective mass in a clean 2D electron system in silicon by analyzing temperature dependence of the weak-field Shubnikov-de Haas (SdH) oscillations in the low-temperature limit rem3 . The effective mass is found to be strongly increased (by a factor of $`3`$) at low electron densities. Using tilted magnetic fields, we find that the value of the effective mass does not depend on the degree of spin polarization, which points to a spin-independent origin of the effective mass enhancement. This is in clear contradiction with existing theories renorm ; spivak ; dol .

Measurements were made in a rotator-equipped Oxford dilution refrigerator with a base temperature of $`30`$ mK on high-mobility (100)-silicon samples similar to those previously used in Ref. krav00 . The resistance, $`R_{xx}`$, was measured by a standard 4-terminal technique at a low frequency (0.4 Hz) to minimize the out-of-phase signal. Excitation current was kept low enough (0.1-0.2 nA) to ensure that measurements were taken in the linear regime of response. Contact resistances in these samples were minimized by using a split-gate technique that allows one to maintain a high electron density in the vicinity of the contacts (about $`1.5\times 10^{12}`$ cm-2) regardless of its value in the main part of the sample. Below, we show results obtained on a sample with a peak mobility close to 3 m2/Vs at $`T=0.1`$ K.

A typical temperature dependence of the amplitude, $`A`$, of the weak-field (sinusoidal) SdH oscillations for the normalized resistance, $`R_{xx}/R_0`$ (where $`R_0`$ is the average resistance), is displayed in Fig. 1. To determine the effective mass, we use the method of Ref. smith  extending it to much lower electron densities and temperatures rem . We fit the data for $`A(T)`$ using the formula

$`A(T)`$ $`=`$ $`A_0{\displaystyle \frac{2\pi ^2k_BT/\mathrm{}\omega _c}{\mathrm{sinh}(2\pi ^2k_BT/\mathrm{}\omega _c)}},`$
$`A_0`$ $`=`$ $`4\mathrm{exp}(2\pi ^2k_BT_D/\mathrm{}\omega _c),`$ (1)
where $`\omega _c=eB_{}/mc`$ is the cyclotron frequency and $`T_D`$ is the Dingle temperature. As the latter is related to the level width through the expression $`T_D=\mathrm{}/2\pi k_B\tau `$ (where $`\tau `$ is the elastic scattering time) afs , damping of the SdH oscillations with temperature may be influenced by temperature-dependent $`\tau `$. We have verified that in the studied low-temperature limit for electron densities down to $`1\times 10^{11}`$ cm-2, possible corrections to the mass value caused by the temperature dependence of $`\tau `$ (and hence $`T_D`$) are within our experimental uncertainty which is estimated at about 10%. Note that the amplitude of the SdH oscillations follows the calculated curve down to the lowest achieved temperatures, which confirms that the electrons were in a good thermal contact with the bath and were not overheated.

In Fig. 2, we show the so-determined effective mass in units of the band mass, $`m_b=0.19m_e`$ (where $`m_e`$ is the free electron mass), as a function of electron density. The effective mass sharply increases with decreasing $`n_s`$. This is in good agreement with the data obtained by an alternative method based on measurements of $`B_c`$ and of the slope of the metallic temperature dependence of conductivity in zero magnetic field, similar to those described in Refs. sh ; sh1 . Using that method, we also determined that our $`g`$ factor is practically density-independent and is equal to $`g2.8`$, which is close to $`g=2`$ in bulk silicon. Note that the agreement between the results obtained using two independent methods as well as the fact that the experimental dependence $`A(T)`$ follows the theoretical curve justify applicability of Eq. (1) to this strongly-interacting electron system.

A strong enhancement of $`m`$ at low $`n_s`$ may originate from spin effects renorm ; spivak ; dol . With the aim of probing a possible contribution from the spin effects, we introduced a parallel magnetic field component, $`B_{}`$, to align the electrons’ spins. As the thickness of the 2D electron system in Si MOSFETs is small compared to the magnetic length in accessible fields, the parallel field couples largely to the electrons’ spins while the orbital effects are suppressed sim ; rem1 .

In Fig. 3, we show the main result of the paper: the behavior of the effective mass with the degree of spin polarization, $`p=(B_{}^2+B_{}^2)^{1/2}/B_c`$. Within our accuracy, the effective mass $`m`$ does not depend on $`p`$. Therefore, the $`m(n_s)`$ dependence is robust, the origin of the mass enhancement has no relation to the electrons’ spins and exchange effects rem2 .

In Fig. 4(a), we compare the extracted Dingle temperature, $`T_D`$, with that recalculated from the electron lifetime determined from zero-field mobility sh1 . Although the elastic scattering time defining the Dingle temperature is in general different from the transport scattering time, the two $`T_D(n_s)`$ dependences are consistent with each other indicating dominant large-angle scattering afs . This indicates that the Dingle temperature decreases with decreasing electron density, in agreement with the narrowing of the cyclotron resonance line observed at low $`n_s`$ in Si MOSFETs cr . In contrast, Figure 4(b) shows that with increasing degree of spin polarization, the Dingle temperature remains approximately constant, while the resistance increases by a factor of about five. So, the straightforward correlation between the Dingle temperature and the elastic scattering time does not hold for fully spin-polarized electron system.

We now discuss the results obtained for the effective mass. We stress that in the present case, the interaction parameter, $`r_s`$, is larger by a factor of $`2m/m_b`$ than the Wigner-Seitz radius and reaches approximately 50, which is above the theoretical estimate for the onset of Wigner crystallization. As has already been mentioned, two approaches to calculate the renormalization of $`m`$ and $`g`$ have been formulated. The first one exploits the Fermi liquid model extending it to relatively large $`r_s`$. Its main outcome is that the renormalization of $`g`$ is large compared to that of $`m`$ renorm . In the limiting case of high $`r_s`$, one may expect a divergence of the $`g`$ factor that corresponds to the Stoner instability. These predictions are in obvious contradiction to our data: (i) the dilute system behavior in the regime of the strongly enhanced susceptibility — close to the onset of spontaneous spin polarization and Wigner crystallization — is governed by the effective mass, rather than the $`g`$ factor, through the interaction parameter $`r_s`$; and (ii) the insensitivity of the effective mass to spin effects also cannot be accounted for.

The other theoretical approach either employs analogy between a strongly interacting 2D electron system and He3 spivak  or applies Gutzwiller’s variational method br  to Si MOSFETs dol . It predicts that near the crystallization point, the renormalization of $`m`$ is dominant compared to that of $`g`$ and that the effective mass diverges at the transition. Although the sharp increase of the mass is in agreement with our findings, it is the expected dependence of $`m`$ on the degree of spin polarization that is not confirmed by our data: the model of Ref. spivak  predicts that the effective mass should increase with increasing spin polarization, whereas the prediction of the other model dol  is the opposite.

Thus, the existing theories fail to explain our finding that in a dilute 2D electron system the effective mass is strongly enhanced and does not depend on the degree of spin polarization. The fact that the spin exchange is not responsible for the observed mass enhancement reduces somewhat the chances for the occurrence of an intermediate phase — ferromagnetic Fermi liquid — that precedes electron crystallization. In principle, should the spin exchange be small, the spin effects may still come into play closer to the onset of Wigner crystallization where the Fermi energy may continue dropping as caused by mass enhancement.

In summary, we have found that in a dilute 2D electron system in silicon, the strongly enhanced effective mass is independent of the degree of spin polarization. This shows that the mechanism underlying the effective mass enhancement is not related to spin and exchange effects.

We gratefully acknowledge discussions with I. L. Aleiner, D. Heiman, and B. Spivak. This work was supported by NSF grant DMR-9988283, the Sloan Foundation, the RFBR, the Russian Ministry of Sciences, and the Programme “The State Support of Leading Scientific Schools”.