text:
# Effects of Al doping on the structural and electronic properties of Mg1-xAlxB2

## Abstract

We have studied the structural and electronic properties of $`\mathrm{Mg}_{1x}\mathrm{Al}_x\mathrm{B}_2`$ within the Virtual Crystal Approximation (VCA) by means of first-principles total-energy calculations. Results for the lattice parameters, the electronic band structure, and the Fermi surface as a function of Al doping for $`0x0.6`$ are presented. The ab initio VCA calculations are in excellent agreement with the experimentally observed change in the lattice parameters of Al doped $`\mathrm{MgB}_2`$. The calculations show that the Fermi surface associated with holes at the boron planes collapses gradually with aluminum doping and vanishes for $`x=0.56`$. In addition, an abrupt topological change in the $`\sigma `$-band Fermi surface was found for $`x=0.3`$. The calculated hole density correlates closely with existing experimental data for $`T_c(x)`$, indicating that the observed loss of superconductivity in $`\mathrm{Mg}_{1x}\mathrm{Al}_x\mathrm{B}_2`$ is a result of hole bands filling.

The discovery of superconductivity in the simple binary intermetallic compound $`\mathrm{MgB}_2`$ with a $`T_c`$ as high as 40 K has stimulated intense investigations, both from the experimental and the theoretical points of view. The superconducting transition temperature $`T_c`$ for $`\mathrm{MgB}_2`$ has been studied as a function of pressure and alloying. Pressure studies have shown that $`T_c`$ decreases with applied hydrostatic pressure, which has been explained by an increase of the band filling of the boron $`\sigma `$-bands with pressure. Thus, the change in the band filling in $`\mathrm{MgB}_2`$ under pressure is an effect of the reduction of the cell volume with pressure.

Experimentally it has been observed that the superconducting transition temperature of $`\mathrm{Mg}_{1x}\mathrm{Al}_x\mathrm{B}_2`$ decrease with Al doping , and superconductivity disappears for $`x>0.5`$ . According to band structure calculations of $`\mathrm{MgB}_2`$, electron doping reduces the density of states (DOS) at the Fermi level. Based on the Rigid Band Approximation, An and Pickett analyzed the effect of Al doping on the DOS of $`\mathrm{MgB}_2`$, and found that the DOS at the Fermi level drops for $`x0.25`$. Structural characterization of $`\mathrm{Mg}_{1x}\mathrm{Al}_x\mathrm{B}_2`$ shows that the cell volume also decreases with Al doping. Therefore, the $`\sigma `$-band filling in $`\mathrm{Mg}_{1x}\mathrm{Al}_x\mathrm{B}_2`$ is expected to have two contributions: the first is due to the electron doping, and the second is a result of the cell volume reduction as in the case of pressure effects.

Measurements of the thermoelectric power, $`S`$, on $`\mathrm{Mg}_{1x}\mathrm{Al}_x\mathrm{B}_2`$ for $`x0.1`$ have shown that the slope of the linear part of $`S(T)`$ changes with Al doping, indicating changes in the Fermi surface due to electron doping. More recently, it was shown that the Raman spectra of $`\mathrm{Mg}_{1x}\mathrm{Al}_x\mathrm{B}_2`$ for $`0x0.6`$ show a pronounced frequency shift and a considerable change in the line-width for the $`E_{2g}`$ phonon mode at $`x0.3`$, which correlate with a steeping in the behavior of $`T_c(x)`$ with Al doping. Although some of the observed effects of Al doping in $`\mathrm{Mg}_{1x}\mathrm{Al}_x\mathrm{B}_2`$ can be interpretated qualitatively in terms of the Rigid Band Approximation as an effect of the $`\sigma `$-band filling, a quantitative analysis is essential in order to determine the interplay between electron doping and the structural, electronic, and transport properties of $`\mathrm{Mg}_{1x}\mathrm{Al}_x\mathrm{B}_2`$.

In this paper, we present a study of the effects of Al doping on the structural and electronic properties of $`\mathrm{Mg}_{1x}\mathrm{Al}_x\mathrm{B}_2`$ for $`0x0.6`$, using the ab initio Virtual Crystal Approximation. The calculated lattice parameters are compared with available experimental data. The evolution of the electronic band structure and the $`\sigma `$-band Fermi surface (FS) as a function of Al doping is analyzed. We correlate our results to the experimentally observed behavior of $`T_c`$ with Al doping. We show that the observed loss of superconductivity in $`\mathrm{Mg}_{1x}\mathrm{Al}_x\mathrm{B}_2`$ can be explained by the filling of the hole bands.

The Khon-Sham total energies were calculated self-consistently using the full-potential linearized augmented plane-wave method (LAPW) as implemented in the WIEN97 code, where the core states are treated fully relativistically, and the semicore and valence states are computed in a scalar relativistic approximation. The exchange correlation potential was evaluated within the generalized gradient approximation (GGA), using the recent parameter-free GGA form by Perdew, Burke, and Ernzerhof. We chose muffin-tin radii $`(R_{MT})`$ of 1.8 and 1.5 a.u. for Mg and B, respectively, and used a plane-wave cutoff $`R_{MT}K_{MAX}=9.0`$. Inside the atomic spheres the potential and charge density were expanded in crystal harmonics up to $`l=10`$. Convergence was assumed when the energy difference between the input and output charge densities was less than $`1\times 10^5`$ Ryd. Special attention was paid to convergence of results by performing the calculations for a sufficiently large number of $`k`$ points in the irreducible wedge of the Brillouin zone for the omega structure (144 $`k`$ points). The corrected tetrahedron method was used for Brillouin-zone integration.

The Al doping was modeled in the ab initio Virtual Crystal Approximation (VCA). The Mg $`(Z=12)`$ sites are substituted by pseudo-atoms which have a fractional electronic charge $`(Z=12+x)`$, depending on the Al concentration, $`x`$. This approximation is justified mainly by the fact that Al only has one electron more than Mg. The full-potential for the VCA system it is determined self-consistently for each value of Al doping without shape approximation. The ab initio VCA as implemented in this work has been used very recently to model C, Cu and Be substitutions in MgB2. The equilibrium lattice parameters were determined by total-energy calculations for each value of Al doping ($`x`$=0.0, 0.1, 0.2, 0.25, 0.3, 0.35, 0.4, 0.5, and 0.6). Since the AlB2 (omega) structure has two structural parameters ($`a`$ and $`c`$), we performed self-consistent total-energy calculations for nine different volumes and for nine different $`c/a`$ ratios, in order to optimize both $`V`$ and $`c/a`$ for each Al concentration. For MgB2 we have obtained $`a=3.083`$ Å and $`c=3.526`$ Å, which compares very well to the experimental values of $`a=3.086`$ Å and $`c=3.524`$ Å , respectively.

In Fig. 1 we present the calculated lattice parameters ($`a`$ and $`c`$) of $`\mathrm{Mg}_{1x}\mathrm{Al}_x\mathrm{B}_2`$ for $`0x0.6`$, and it can be seen that both $`a`$ and $`c`$ decrease with Al doping as observed experimentally. For comparison we have included the experimental data from Slusky et al., and we find that the change in $`a`$ with Al doping is very well reproduced by the VCA calculations. For the region $`0.1x0.2`$ two values of $`c`$ were reported, which has been ascribed to the coexistence of two phases (Mg-rich and Al-rich phases). It is interesting to note that in this region the VCA values approximately reproduce the average value. However, for $`x>0.2`$ we find very good agreement between the experimental data and the VCA calculations (see Fig. 1). Although both cell parameters ($`a`$ and $`c`$) decrease monotonically with Al doping, it is interesting to note that the slope for $`c`$ as a function of Al doping is larger than for $`a`$. In order to understand this behavior in the context of bonding properties we have analyzed the change in the charge distribution with Al doping. Fig. 2(a) shows the charge density distribution in the (110) plane of MgB2. Mg nuclei are located at the corners of the map and B nuclei are at the $`(1/3,1/2)`$ and $`(2/3,1/2)`$ positions, all of them in the plane of the figure. We can see the directional, covalent B-B $`\sigma `$-bonds. In addition, there is a significant density of charge in the interstitial region giving rise to metallic-type bonding between the Mg and B planes. The charge distribution and bonding properties of MgB2 have been calculated previously and were discussed in detail in Ref. 19. Therefore, we concentrate on the influence of Al doping on the electron density of $`\mathrm{MgB}_2`$. Fig. 2(b) shows the difference between the charge densities of $`\mathrm{Mg}_{0.5}\mathrm{Al}_{0.5}\mathrm{B}_2`$ and $`\mathrm{MgB}_2`$. It can be seen that charge transfer occurs from Al-B ions into the nearby interstitial region. We can see that the majority of this charge is distributed in the inter-plane region, and an important fraction of the charge is being transfered to the $`\pi `$-bond, while only a small fraction is at the $`\sigma `$-bond in the boron planes. This important increase in occupation of the $`\pi `$-bond with Al doping accounts for the strong decrease of the separation between planes (the $`c`$ axis) of $`\mathrm{Mg}_{1x}\mathrm{Al}_x\mathrm{B}_2`$.

The evolution of the calculated DOS for $`\mathrm{Mg}_{1x}\mathrm{Al}_x\mathrm{B}_2`$, not shown here, shows a band broadening as a function of Al doping, mainly as a consequence of the cell volume reduction. In addition, electron doping raises the Fermi level to higher energies. Both effects, band broadening and electron doping, contribute to reduce the density of states at the Fermi level in Al-doped $`\mathrm{MgB}_2`$. In Table I, we summarize the calculated lattice parameters, cell volume, and the total density of states at the Fermi level \[$`N(E_F)`$\] for each of the studied Al concentrations. We can see that $`N(E_F)`$ decreases with Al doping, from 0.72 for $`\mathrm{MgB}_2`$ to 0.26 states/eV per cell for $`\mathrm{Mg}_{0.4}\mathrm{Al}_{0.6}\mathrm{B}_2`$. Therefore, in a BCS scenario this reduction in $`N(E_F)`$ accounts for the decrease of $`T_c`$ with Al doping in $`\mathrm{Mg}_{1x}\mathrm{Al}_x\mathrm{B}_2`$.

A careful analysis of the $`x`$-dependence of the electronic band structure, and in particular of the $`\sigma `$-band FS which has been shown to be relevant for superconductivity in $`\mathrm{MgB}_2`$, provides a more detailed and quantitative description of the effects of Al doping. In Fig. 3 we present the electronic band structure of $`\mathrm{MgB}_2`$. The $`\sigma `$-bands coming from the $`s`$-$`p`$ boron orbitals, are strongly two-dimensional with very little dispersion along $`\mathrm{\Gamma }`$-$`A`$, this dispersion can be characterized by the difference between the $`E_\mathrm{\Gamma }`$ and $`E_A`$ energies (see Fig. 3). The $`E_\mathrm{\Gamma }`$ and $`E_A`$ energies correspond to the bottom and top of the $`\sigma `$-band in the $`\mathrm{\Gamma }`$-$`A`$ direction, respectively. The light-hole and heavy-hole $`\sigma `$-bands in $`\mathrm{MgB}_2`$ form a FS consisting of two fluted cylinders surrounding the $`\mathrm{\Gamma }`$-$`A`$ line in the Brillouin zone (see Fig. 3). The dependence of the energy of the $`\sigma `$-bands at $`\mathrm{\Gamma }`$ and $`A`$ relative to $`E_F`$ as a function of Al doping are shown in Fig. 4. We can see that both energies, $`E_\mathrm{\Gamma }`$ and $`E_A`$, decrease monotonically as a function of Al doping. More interestingly, the Fermi level reaches $`E_\mathrm{\Gamma }`$ for $`x=0.3`$ and $`E_A`$ for $`x=0.56`$. We find that the radius of the cylinders decreases gradually with Al doping and at a critical concentration of $`x=0.3`$, the radius at $`k_z=0`$ collapses and the FS takes the form of a sandglass. A three-dimensional view of the changes in the FS topology with Al doping are presented in Fig. 4. For $`x=0.3`$ the Fermi level in $`\mathrm{\Gamma }`$ is at a saddle point in the band structure, and the transition through the saddle point results in the disruption of the neck, i.e., the transition from a closed to an open section of FS. For $`x>0.3`$ the FS takes the form of two cones (see Fig. 4), and these finally vanishes at the second critical concentration $`(x=0.56)`$ when the hole bands have been filled. These changes in the hole FS are expected to be accompanied by various kinds of electronically driven anomalies, including lattice dynamics and transport properties.

As was mentioned above, Raman spectroscopy studies on $`\mathrm{Mg}_{1x}\mathrm{Al}_x\mathrm{B}_2`$  show a pronounced shift and a considerable change in the line-width of the $`E_{2g}`$ phonon mode at $`x0.3`$. Additionally, a steeping of the $`T_c`$ decrease has been observed at an Al concentration of approximately 0.3. These changes in both the structural and the superconducting properties seems to be strongly related to the abrupt change in the FS topology, which occurs for $`x=0.3`$ (see Fig. 4). In order to establish a more direct comparison between the FS evolution and the superconducting properties with increasing the Al doping, we have calculated the hole FS area as a function of $`x`$, which is proportional to the hole density at the Fermi level. In Fig. 4(b) we compare the calculated normalized FS area, $`A_{FS}(x)/A_{FS}(0)`$, with the normalized superconducting critical temperaure, $`T_c(x)/T_c(0)`$. The experimental data for $`T_c(x)`$ were taken from Ref. 8. We can see that for the low concentration region ($`x0.25`$), before the $`E_{2g}`$ phonon frequency shift, the drop of $`T_c`$ is directly related to the change in the hole density. This view is in agreement with recent results of NMR experiments on Al-doped $`\mathrm{MgB}_2`$ for $`x0.1`$. In the high concentration region $`(x>0.25)`$, the behavior of $`T_c`$ is determined by the FS area but the importance of the phonon-renormalization is clear. In this way, the FS area and $`T_c`$ follow the same behavior with Al doping in the whole range $`(0x0.6)`$, indicating a close relation between the changes in the $`\sigma `$-band FS and the loss of superconductivity in Al doped $`\mathrm{MgB}_2`$.

In summary, we have performed a first-principles study of the effects of Al doping on the structural parameters, the electronic structure, and the $`\sigma `$-band FS of $`\mathrm{Mg}_{1x}\mathrm{Al}_x\mathrm{B}_2`$, using the Virtual Crystal Approximation. (i) We find that the ab initio VCA calculations are in excellent agreement with the experimentally observed changes in the lattice parameters as a function of Al doping. (ii) The analysis of the charge density shows that an important portion of the Al-electrons are at the inter-plane region and only a small fraction at the B-B planes, providing an explanation for the strong change of the $`c`$-axis and the small change in the $`a`$-axis with the Al concentration. (iii) The hole FS gradually collapses with Al doping and vanishes for $`x=0.56`$. An abrupt topological change was found for $`x=0.3`$, which correlates with the frequency shift of the $`E_{2g}`$ phonon mode and the steeping in the $`T_c(x)`$ decrease. Additionally, the critical concentration of $`x=0.56`$ at which the hole FS disappears, corresponds to the experimentally observed Al concentration (0.5-0.6) for which $`T_c(x)`$ vanishes. (iv) We find that the behavior of the calculated $`\sigma `$-band FS area with Al doping correlates with the superconducting critical temperature $`T_c(x)`$. Consequently, the observed loss of superconductivity in $`\mathrm{Mg}_{1x}\mathrm{Al}_x\mathrm{B}_2`$, can be explained as a result of the filling of the hole bands.

This research was funded by the Consejo Nacional de Ciencia y Tecnología (CONACYT, México) under Grant No. 34501-E. Two of the authors (O.P. and A.A.) gratefully acknowledge a student fellowship from CONACYT-México. The authors would like to thank Dimitris Papaconstantopoulos and David Singh for valuable discussions.


format: text

text:
# On the magnetotransport of 3D systems in quantizing magnetic field

## I General formalism

The main goal of the present paper concerns the peculiar thermodynamic approach\[11-13\] regarding 3D electron transport in strong magnetic fields. We will follow the argumentation first put forward by Kirby and Laubitz , and then modified in Refs. for magnetotransport in 2D electron(hole) systems. This approach allows one to account for both the Schubnikov-de Haas Oscillations and Integer Quantum Hall Effect modes. Let us consider the isotropic sample subjected in strong magnetic field $`B=B_z`$. Neglecting spin splitting, the 3D energy spectrum yields

$$\epsilon _n=\mathrm{}\omega _c(n+1/2)+\epsilon _{},$$
(1)
where $`n=0,1`$… is the LL number, $`\omega _c=eB/mc`$ is the cyclotron frequency, $`m`$ is the electron effective mass. Then, $`p_z`$ and $`\epsilon _{}=\frac{p_z^2}{2m}`$ are the momentum and the kinetic energy respectively of an electron moving along the magnetic field. The electron motion in x-y plane is quantized, thus results in discrete Landau level set. For simplicity, we further neglect the broadening of the Landau levels. The density of states associated with the certain LL obeys the textbook result $`\mathrm{\Gamma }=1/2\pi l_B^2`$, where $`l_B=(\mathrm{}c/eB)^{1/2}`$ is the magnetic length. In the present paper we mostly interest in the strong quantum limit case, when $`\mathrm{}\omega _c>>kT`$.

In general, the macroscopic current, $`𝐣`$ , and the energy flux, $`𝐪`$ , densities are given by

$`𝐣=\widehat{\sigma }(𝐄\alpha T),`$ (2)
$`𝐪=\left(\alpha T\zeta /e\right)𝐣\widehat{\kappa }T.`$
Here, $`𝐄=\zeta /e`$ is the electric field, $`\mu `$ and $`\zeta =\mu e\phi `$ are the chemical and electrochemical potential respectively, $`\alpha `$ is the thermopower. Then, $`\widehat{\sigma }`$ and $`\widehat{\kappa }=LT\widehat{\sigma }`$ are the conductivity and thermal conductivity tensors respectively, $`L=\frac{\pi ^2k^2}{3e^2}`$ is the Lorentz number. It will be recalled that Eq.(2) is valid for a confined-topology sample for which the diamagnetic surface currents  are taken into account. Both the Einstein and Onsager relationships are satisfied. The movement of the electron along the magnetic field is not quantized, thus results in Drude longitudinal component $`\sigma _{zz}=Ne^2\tau /m`$ of the conductivity tensor, where $`N`$ is the 3D electron density. In what follows we assume the constant momentum relaxation time $`\tau `$ in z-direction. Compared to unperturbed movement of an electron along the magnetic field, the transverse drift of electrons in x-y plane results in off-diagonal Hall components $`\sigma _{yx}=\sigma _{xy}=Nec/B`$ of the conductivity tensor. Then, in contrast to previous studies\[2-5\] we suggest that in strong quantum limit $`\sigma _{xx}=\sigma _{yy}=0`$. Actually, the x-y plane related fraction of Eq.(2) can be viewed as the current and energy dissipationless fluxes caused by electron drift in crossed fields. With the help of the above notations, the resistivity tensor $`\widehat{\rho }=\widehat{\sigma }^1`$ obeys the same symmetry, namely $`\rho _{yx}=\rho _{xy}=1/\sigma _{yx}`$, $`\rho _{zz}=1/\sigma _{zz}`$. The crucial idea of the present paper is that a nonzero transverse magnetoresistivity $`\rho `$ can, nethereless, arise due to combined action of the Peltier and Seebeck thermoelectric effects\[11-13\].

## II Results and discussion

Let us discuss in details the typical experimental setup allowed to measure the Hall resistivity $`\rho _{yx}`$ and transverse magnetoresistivity $`\rho `$. The standard Hall-bar geometry sample (Fig.1, inset) is connected to the current source by means of two identical leads. The contacts ”a” and ”b” are ohmic. The longitudinal voltage is measured between the open ends (”e” and ”d”), maintained at the ambient temperature. The sample is placed in a chamber ( not shown in Fig.1, inset) kept at the bath temperature $`T_0`$. Note that the electron temperature may, in principle, exceed the bath temperature if the electron-phonon coupling is weak at low temperatures. Actually, the 3D electron gas cooling can occur predominantly through the contacts of the sample and the leads connected to them. However, the heat leakage via metal leads can be disregarded similar to that reported for two-dimensional electron case. Finally, we will consider the adiabatic cooling of 3D electrons.

Recall that the Peltier heat is generated by a current across the interface between two different conductors. At the contact ( e.g. ”a” in Fig.1, inset), the temperature, $`T_a`$, electrochemical potential $`\zeta `$, normal components of the total current, $`I`$, and the total energy flux are continuous. There exists a difference $`\mathrm{\Delta }\alpha =\alpha _m\alpha `$ between the thermoelectric powers of the metal conductor and 3D sample, respectively. Then, $`Q_a=I\mathrm{\Delta }\alpha T_a`$ is the amount of Peltier heat released per unit time in contact ”a”. For $`\mathrm{\Delta }\alpha >0`$ and current flow direction shown in Fig.1, the contact ”a” is heated and contact ”b” is cooled. The contacts are at different temperatures, and $`\mathrm{\Delta }T=T_aT_b>0`$. At small currents, the temperature gradient is small and $`T_{a,b}T_0`$. In this case the 3D sample thermopower $`\alpha `$ is nearly constant, hence, one can disregard the Thompson heating $`IT\alpha `$ in 3D bulk. Note that the amount of the Peltier heat evolved at the contact ”a” is equal to that absorbed at the contact ”b”. The energy flux is continuous at each contact, thus Eq.(2) yields $`\kappa _{yx}_xT|_{a,b}=j\mathrm{\Delta }\alpha T_{a,b}`$. Here, we take into account that the current is known to enter and leave the sample at two diagonally opposite corners ( Fig.1, inset). One can fulminantly find the longitudinal temperature gradient $`_xT=\frac{j\mathrm{\Delta }\alpha }{L\sigma _{yx}}`$, which is linear in current. Omitting the contribution of the conductor resistances, the voltage drop, $`U`$, measured between the ends ”e” and ”d” is equal to Seebeck thermoelectromotive force $`U=E_x𝑑x=\alpha 𝑑T=\mathrm{\Delta }\alpha (T_aT_b)`$, thus gives the transverse magnetoresistivity $`\rho =U/jl=(\mathrm{\Delta }\alpha )^2\rho _{yx}/L`$.

We now intend to find the actually measured longitudinal resistivity taking into account the Peltier-Seebeck thermoelectric effects. Since the longitudinal movement of an electron unaffected by magnetic field, we make use of the result valid for $`B=0`$. The longitudinal resistivity of 3D electron gas acquire the correction $`\mathrm{\Delta }\rho _{zz}=(\mathrm{\Delta }\alpha )^2\rho _{zz}/L`$. Finally, we summarize the actually measured values of the 3D electron gas resistivities as it follows

$$\rho _{yx}=\frac{B}{Nec},\rho =s\rho _{yx},\rho _{zz}=(1+s)\frac{\rho _{yx}}{\omega _c\tau },$$
(3)
where we take into account that for the actual case of the metal leads $`s=(\mathrm{\Delta }\alpha )^2/L\alpha ^2/L`$. Remind that in the strong quantum limit, the thermopower of dissipationless 3D electron gas is a universal thermodynamic quantity proportional to the entropy per electron:

$$\alpha =\frac{S}{eN},$$
(4)
where $`S=\left(\frac{\mathrm{\Omega }}{T}\right)_{\mu ,B}`$ is the 3D electron entropy, $`N=\left(\frac{\mathrm{\Omega }}{\mu }\right)_{T,B}`$ is the 3D electron density, $`\mathrm{\Omega }=2kT\mathrm{\Gamma }\underset{n,p_z}{}\mathrm{ln}(1+\mathrm{exp}((\mu \epsilon _n)/kT))`$ is the thermodynamic potential, which accounts the spin degeneracy.

Ii is instructive to discuss the applicability of the Gibbs statistics formalism for actual 3D electrons case. The make use of thermodynamic potential presumes a variable number of 3D particles while the chemical potential is assumed to be a constant. We argue that this method is not prohibited since the chemical potential of 3D electrons can be, for example, pined by an external reservoir of carriers. The metal leads can play the role of the reservoir. Indeed, in thermodynamic equilibrium the chemical potential of the multiple system ”3D electrons-conductor” flattens out. Let us suppose a certain value of the chemical potential $`\mu `$, and, thus 3D density $`N_0`$ at $`B=0`$. The magnetic field growth may result in subsequent changes of 3D electron density $`N`$. In strong fields the only few Landau levels are filled, hence, the absolute change in $`N`$ can be of the order of the zero-field carrier density $`N_0`$. Simultaneously, the number electrons in the external reservoir changes. The relative variation of the chemical potential of the external reservoir yields $`\delta \mu /\mu =\frac{2}{3}N_0/N_m10^7`$, where $`N_m10^{23}`$cm-3 is the electron density the metal leads, $`N_0=10^{16}`$cm-3 is the typical density of 3D solid. We conclude that the chemical potential of metal leads, and, hence 3D electrons remains constant. Note that conventional reasoning of 3D electron plasma neutrality are rather doubtful since the Debye screening length can be of the order of macroscopic sample length scale( see Appendix IV ) in quantizing magnetic field.

It can be demonstrated that both the electron density $`N`$ and the entropy $`S`$ are universal functions of the reduced temperature $`\xi =kT/\mu `$ and dimensionless magnetic field $`\mathrm{}\omega _c/\mu =\nu ^1`$, where $`\nu `$ is the so-called filling factor. Following the conventional Lifshitz-Kosevich formalism, we can easily derive asymptotic equations for $`N,S`$, and, hence, for $`\alpha `$, valid in the low-temperature and low-field limit $`\nu ^1,\xi <1`$:

$`N`$ $`=`$ $`N_0\left(1+{\displaystyle \frac{\pi ^2}{2F_{1/2}\left(\frac{1}{\xi }\right)}}{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^k\mathrm{sin}\left(2\pi k\nu \frac{\pi }{4}\right)}{\sqrt{z}\mathrm{sinh}z}}\right),`$ (5)
$`S`$ $`=`$ $`S_0+{\displaystyle \frac{N_0k\pi ^3}{2F_{1/2}\left(\frac{1}{\xi }\right)}}{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}(1)^k\mathrm{\Phi }(z)\mathrm{cos}\left(2\pi k\nu {\displaystyle \frac{\pi }{4}}\right),`$
where $`\mathrm{\Phi }(z)=\frac{z\mathrm{coth}z1}{z^{3/2}\mathrm{sinh}z}`$ is a form-factor, $`z=2\pi ^2\xi \nu kkT/\mathrm{}\omega _c`$ is the dimensionless parameter which roughly corresponds to quantum limit criteria at $`k=1`$. Then, $`N_0=n_0\frac{3}{2}\xi ^{3/2}F_{1/2}(1/\xi )`$ and $`S_0=N_0k\left(\frac{5F_{3/2}(1/\xi )}{3F_{1/2}(1/\xi )}\frac{1}{\xi }\right)`$ are the zero-field electron gas density and entropy respectively, $`F_n(y)`$ is the Fermi integral and $`n_0=\frac{(2m\mu )^{3/2}}{3\pi ^2\mathrm{}^3}`$ is the density of strongly degenerate 3D electron gas at $`T=0`$.

In Fig.1 we plot the magnetic field dependencies of the Hall resistivity $`\rho _{yx}(\nu ^1)`$ and transverse magnetoresistivity $`\rho (\nu ^1)`$ given by Eqs.(3-5). Both curves are scaled in units of $`\rho _c=\frac{m\mu }{n_0e^2\mathrm{}}=\frac{h}{e^2}\frac{3\pi }{4k_F}`$, where $`k_F=\sqrt{2m\mu }/\mathrm{}`$ is the Fermi vector. Then, the longitudinal resistivity $`\rho _{zz}(\nu ^1)`$ specified by Eq.(3) is scaled in units of zero-field resistivity $`\rho _{zz}(0)=\frac{m}{n_0e^2\tau }=\frac{2\rho _c}{k_Fl}`$, where $`l=`$ $`\mathrm{}k_F\tau /m`$ is the mean free path. The overall behavior of the resistivity curves shown in Fig.1 can be easily understood in terms of the 3D energy spectrum, and, thus the related density of states. Indeed, apart from the discrete Landau levels, the 3D energy spectrum specified by Eq.(1) contains the continuous component $`\epsilon _{}`$. At fixed $`\mu `$ the number of occupied states, and, hence $`\rho _{yx},\rho _{zz}N^1`$ are smooth functions of the magnetic field. Then, the density of states is continuous, hence provides the nonzero transverse magnetoresistivity $`\rho `$. It is worthwhile to mention that in 2D electron gas case the energy spectrum is purely discrete, thus leads to quantized Hall resistivity $`\rho _{yx}`$, with the longitudinal resistivity $`\rho `$ vanishing within the Hall plateaus.

We now make an attempt to compare our results with experimental observations. As expected, the longitudinal, transverse and the Hall components of magnetoresistivity exhibit familiar 1/H-period in-phase oscillations associated with discrete LL energy spectrum. However, in contrast to our predictions the behavior of the longitudinal and the transverse magnetoresistivity components becomes uncorrelated at high fields when the only few LLs are filled. The possible reason for the above discrepancy could be, for example, the constant momentum relaxation time, i.e. $`\tau \tau (B)`$, adopted in our simple model. For typical n-InSb sample $`n_0=1.210^{16}`$cm-3 used in Ref. we obtain $`\mu =170`$K and, then $`\rho _c=0.1`$Om$``$cm. At filling $`\nu =3/2`$ our zero spin-split approach provides(see Fig.1) the transverse magnetoresistivity amplitude $`\rho 0.002`$ Om$``$cm at liquid helium temperature $`\xi =0.025`$. This value is consistent with that $`0.015`$ Om$``$cm observed experimentally at $`B=2.3`$T.

## III Conclusions

In conclusion, we have calculated the longitudinal, transverse and the Hall magnetoresistivities of a 3D solid placed in quantizing magnetic taking into account the contribution caused by combined action of Peltier and Seebeck thermoelectric effects. The transport coefficients demonstrate familiar 1/H-period oscillations being universal functions of the filling factor. The amplitude of the transverse resistivity is consistent with that observed in experiment.

## IV Appendix

Let us examine the charge relaxation and screening for 3D solids placed in magnetic field. These effects can be described by continuity equation

$$\frac{\rho _e}{t}+div𝐣=0.$$
(6)
In quantizing magnetic field the non-dissipative transverse current specified by Eq.(2) can be re-written as it follows $`𝐣=\frac{cN}{B^2}`$rot$`(\zeta 𝐁)`$. Thus, the charge relaxation is absent since div$`𝐣0`$, therefore $`\rho _e/t=0`$. The retardation effects, if accounted, could result in slow charge relaxation $`(\sigma _{yx}/c)^2`$ analogous to that in 2D electron gas case. We argue that the transverse length of the Debye screening could be comparable with the macroscopic( i.e. sample) length scale. We conclude that in quantizing magnetic field the electron plasma neutrality may be violated in 3D sample bulk.


text:
# Observation of a new excitation in bcc 4He by inelastic neutron scattering

## Abstract

We report neutron scattering measurements of the phonons in bcc solid 4He. In general, only 3 accoustic phonon branches should exist in a monoatomic cubic crystal. In addition to these phonon branches, we found a new ”optic-like” mode along the  direction. One possible interpretation of this new mode is in terms of localized excitations unique to a quantum solid.

The main attribute of a quantum solid is the large zero point motion of the atoms. The vibrating atoms thus encounter the repulsive part of the interatomic potential, leading to strong short range correlations. Theoretical description of the solid is consequently a tour de force study of anharmonic effects. The self consistent phonon (SCP) theory which emerged over the years, treats these short range correlations using Jastrow-type wavefunctions in the description of the ground state, in addition to the inclusion of anharmonic terms in the potential. Phonon dispersion curves calculated in the SCP theory are in good agreement with neutron scattering experiments on solid 4He , except for the T1 branch of bcc 4He, which experimentally has a much lower energy. In addition, the SCP theory was able to explain several observations unique to quantum solids, such as multiphonon scattering which dominates the scattering already at intermediate momentum transfer Q, and phonon interference which gives rise to oscillations of the effective Debye Waller factor. Our motivation to do new experiments came following a recently proposed model . This model aims to explain the discrepancy regarding the T1 branch of bcc 4He through an approach in which the zero point motion is treated separately from the usual harmonic degrees of freedom. Indeed, the T1 mode predicted by the new model agrees with experiment, while leaving the other phonons unaffected. In addition, the model predicts the appearance of a new dispersionless excitation branch at an energy twice that of the T1 phonon at the zone boundary ($``$1.2meV). Such an ”optic-like” mode was not reported.

In order to check this prediction, we performed inelastic neutron scattering experiments on the bcc phase of 4He using the IN-12 spectrometer at the ILL. 4He single crystals were grown and oriented in the beam. Since the maximal temperature width of the bcc phase is less than 50 mK, particular care was devoted to temperature control, resulting in a $`\pm `$0.1mK stability. In addition, the growth cell contained a capacitive pressure gauge allowing us to precisely control the location on the phase diagram. Out of the 5 crystals grown, 4 could be oriented in the {11$`\overline{1}`$} scattering plane, and one in the {100}. From the observed Bragg peaks, the cell contained only one large single crystal. The energy of the incident beam could be varied in the range of 2.3 to 14 meV. Constant-Q scans were usually done at fixed kF. In the high resolution scans we used a cooled Be filter to remove $`\lambda /2`$ contamination. The highest instrumental resolution used was 0.07 meV (FWHM) with kF=1.15 1/Å. Data were typically taken around the (1,1,0) point.

In addition to well defined single phonon groups and broad multiphonon features, the measured spectra also show a new, well defined feature. In Fig. 1, we show an intensity map displaying longitudinal scans in which the new feature appears as a weakly dispersive ”optic-like” branch, together with the L phonons. The resolution in Fig. 1 is not sufficient to display details of the spectra, which are seen more clearly in Fig. 2. The new feature was observed in all 4 crystals during scans along the  direction. This feature disappeared once the crystal was molten and the cell contained only liquid. It is therefore not connected with scattering from the liquid or from the walls of the cell . Moreover, we performed intentional experiments in the mixed phase region, with both hcp and bcc crystals in the cell. We observed no correlation between the phonon peaks of the hcp phase and the new feature. We conclude therefore that the new feature is a property of the bcc solid.

The average energy (over all 4 crystals) of this new mode at the zone origin is 1.23 meV, almost exactly twice the energy of the zone boundary T1 phonon (0.59$`\pm `$0.01 meV). This ratio is constant over the whole temperature range where the bcc phase exists. The best location to observe this mode was at a reduced wavevector $`q`$=0, where it was not masked by any of the other phonon branches. However, as Fig. 1 shows, it is possible to resolve this mode up to $`q`$=0.2 in units of $`2\pi /a`$ (r.l.u.). At higher $`q`$, the phonon groups are too broad to identify this feature. In Fig.3, we show the dispersion relations measured along the  direction. The lines in the figure are from the SCP theory  and from the new model . Our data for the three phonon branches is in excellent agreement with previous experiments. The fact that the T1 branch is lower by a factor of $``$ 2 compared with the SCP calculation was already evident from the previous experiments , and is confirmed in the present work. In addition, Fig. 3 shows the data for the new feature, which shows some dispersion. It is important to ask whether the new feature is a distinct mode, rather than a product of anharmonic effects described by the SCP theory, namely phonon interference and multiphonon scattering. First, phonon interference gives rise to scattering intensity which is an odd function of $`q`$, namely it vanishes at the zone origin. The fact that we see this feature clearly at the zone origin excludes the possibility that it is due to interference effects. Regarding multiphonon effects, one could generate a feature at $`q`$=0 and energy of 1.2 meV as a bound state of two zone boundary T1 phonons with opposite momenta. However, in the same way one could combine different phonons to obtain any energy and $`q`$, namely this process will produce a broad background, in contrast to what is seen in Fig  2. Typically, the width of multiphonon features is on the order of $``$several meV, while that of the new feature is an order of magnitude smaller (about 0.2 meV). Next, in Fig. 4 we plot the measured linewidth of all the phonons along (110) vs. their energy. To obtain the data in this figure, we used well defined, symmetric phonon groups. The phonon groups were fit to a damped harmonic oscillator function folded with the gaussian spectrometer resolution function determined by scattering from Vanadium sample. One can see that phonons with energy less than 1.2 meV have widths which are below our resolution limit ($``$0.02 meV), while above that energy there is a marked increase of the linewidth. In contrast, the linewidth predicted by the SCP theory due to phonon interaction is a smoothly increasing function of energy (see Fig. 4). The abrupt increase shown by the data strongly suggests that the feature at 1.23 meV is a distinct excitation which the high energy phonons can decay into.

We outline here possible interpretation of the data involving localized excitations. Here, localized excitations can be point defects or ”local modes”. This approach is appealing since point defects in solid He are predicted to add another excitation branch to the usual phonon spectrum. Interpretation of the new feature in these terms is supported by thermally activated mass  and spin diffusion experiments in bcc 4He, in which the activation energy of point defect was found to be 1.25 meV. One possibility is that the new observed mode is connected to a point defect such as vacancy, interstitial, or a Frenkel pair, involving a displacement of atoms from their lattice sites. Several authors proposed that point defects in quantum crystals are delocalized objects and can behave as quasiparticles, having a finite bandwidth and well defined dispersion relation. Indeed, the data for the new feature in Fig. 2 can be fit by an expression of the type $`E(q)=ϵ_0+\mathrm{}^2q^2/2m^{}`$, with $`ϵ_0`$=1.23 meV and $`m^{}`$ 0.2 m4. We cannot determine the bandwidth directly because the new mode cannot be resolved for $`q`$0.2 r.l.u. There are however problems associated with this interpretation: First, according to the data shown in Fig. 4, phonons can decay into this mode once their energy exceeds $`ϵ_0`$. For the L branch, the phonon with the lowest energy(=$`ϵ_0`$) that can decay has $`q`$=0.15 r.l.u. If the mode has intrinsic dispersion, then there is a problem with momentum conservation. Since the decay is observed, it is possible that the new mode is intrinsically dispersionless, and the apparent dispersion seen at finite $`q`$ can be attributed to mode coupling effects. A dispersionless mode is not what these models predict. Second, these models assume that the energy spectrum of point defects exists along every direction in the crystal. Our preliminary measurements done along the  and  directions do not show the new mode.

The other possibility is an interpretation according to the new model. In this model, this new excitation represents a density fluctuation localized over two unit cells. Unlike the point defects described above, this excitation is a dynamic density fluctuation, in which atoms remain associated with their lattice sites. First, the model predicts correctly the energy of this new excitation at $`q`$=0. Second, as fig. 3 shows, both the L and T2 phonons become damped above the same energy of 1.2 meV. For the L phonon this happens at $`q`$=0.15 r.l.u., while for the T2 phonon the same thing happens at $`q`$=0.27 r.l.u. Hence, the new feature seems to be essentially dispersionless, as predicted. However, one problem with this interpretation is that coupling of the new mode with other phonons, leading to the apparent dispersion and phonon decay, is not treated within the model.

Finally, model calculations done on strongly anharmonic solids indicate the possible existence of localized modes. It would be interesting to extend these calculations using a potential more specific to solid 4He, to allow a comparison with experiment.

In conclusion, we have observed a new excitation in bcc 4He. One interpretation which we suggest views this excitation as a localized excitation unique to a quantum crystal. The nature of the localized excitation relies on the specific way in which atomic zero point motion is incorporated into the theoretical models.

Acknowledgments:  We thank S. Hoida and M. Ayalon (Technion), S. Raymond, F. Thomas, S. Pujol and J. Previtali (ILL) for their help with the experiment. We have benefitted from discussions with B. Fåk, H. R. Glyde, A. Auerbach, B. Shapiro, and S. G. Lipson. This work was supported in part by The Israel Science Foundation, and by the Technion Fund for Promotion of Research.

Figure Captions

Fig. 1 Intensity map of longitudinal scans along  showing the L branch and an ”optic like” excitation band. Note that the intensity scale in the figure is logarithmic, given in neutron counts.

Fig. 2 Typical neutron scattering spectra along  showing the new feature for two values of $`q`$. In (a), taken at $`q`$=0.1, the feature appears together with the L phonon. In (b), taken at the zone origin ($`q`$=0), only the new feature is seen. The lines in the figure are gaussian, and serve as guide to the eye.

Fig. 3 Dispersion curves along . Solid symbols - data for the 3 phonon branches. Open symbols-the new feature. Solid lines represent the SCP theory . The additional predictions of the new model are shown as dash-dotted line (The prediction of the new model for the L and T2 branches are the same as of the SCP theory).

Fig. 4 Experimental phonon linewidth vs. phonon energy for all 3 branches along  showing the abrupt increase near 1.2 meV. Symbols are experimental data. Dash-dot line - SCP theory for the L phonon (ref.). The dashed line is a guide to the eye.


text:
# On a probable manifestation of Hubble expansion at the local scales, as inferred from LLR data

## 1 Introduction

Precise measurements of the Earth–Moon distance by using the ultrashort laser pulses – the lunar laser ranging (LLR) – are carried out for over 30 years, after the installation of several retroreflectors on the lunar surface in the course of Apollo (USA) and Lunakhod (USSR, in collaboration with France) space missions (e.g. reviews by Dickey et al. 1994; Nordtvedt 1999; Samain et al. 1998). The typical accuracy of these measurements was:

* $`25`$ cm in the early 1970’s,
* 2–3 cm in the late 1980’s, and
* a few millimeters at the present time.

LLR works contributed significantly to astrometry, geodesy, geophysics, lunar planetology, and gravitational physics. The most important results related to General Relativity are:

* verification of the Strong Equivalence Principle with accuracy up to $`10^{13}`$,
* determination of the first post-Newtonian parameters in the gravitational field equations,
* detection of the so-called geodetic precession of the lunar orbit, and
* imposing the observational constraints on time variations in the gravitational constant $`\dot{G}/G`$ with accuracy $`10^{11}`$ per year.

## 2 Using the LLR data for finding the “local” Hubble constant

Despite the considerable advances listed above, there is a long-standing unresolved problem in the interpretation of LLR data – anomalous increase in the lunar semimajor axis (e.g. Pertsev 2000).

In general, such increase is well-known and can be *partially* explained by tidal interaction between the Earth and Moon (e.g. Kaula 1968). Because of the relaxation processes, the tidal bulge is not perfectly symmetric about the Earth–Moon line but slightly shifted in the direction of Earth’s rotation (see Fig. 1). As a result, there is a torque moment, which decelerates a proper rotation of the Earth and accelerates an orbital rotation of the Moon; so that the mean Earth–Moon distance increases.

From the angular momentum conservation law

$$I_\mathrm{E}\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{\Omega }_\mathrm{E}+m_\mathrm{M}\frac{\mathrm{d}}{\mathrm{d}t}\left(R^2\mathrm{\Omega }_{\mathrm{ME}}\right)=0$$
(1)
and the relation between the lunar orbital velocity and its distance from the Earth

$$\mathrm{\Omega }_{\mathrm{ME}}=G^{1/2}m_\mathrm{E}^{1/2}R^{3/2},$$
(2)
it can be easily found that the rate of increase in the lunar semimajor axis $`\dot{R}`$ is related to the rate of deceleration of the Earth’s rotation $`\dot{T}_\mathrm{E}`$ by the well-known formula:

$`\dot{R}`$ $`=`$ $`k\dot{T}_\mathrm{E},`$ (3)
$`k`$ $`=`$ $`4\pi G^{1/2}I_\mathrm{E}m_\mathrm{E}^{1/2}m_\mathrm{M}^1R^{1/2}T_\mathrm{E}^2`$ (4)
$`=`$ $`1.8110^5\mathrm{cm}/\mathrm{s},`$
where $`\mathrm{\Omega }_\mathrm{E}`$ and $`T_\mathrm{E}`$ are the angular velocity and period of the proper rotation of the Earth, $`\mathrm{\Omega }_{\mathrm{ME}}`$ is the angular velocity of orbital rotation of the Moon about the Earth, $`R`$ is the distance between them111Within the accuracy required here, we can neglect the ellipticity of the lunar orbit., $`I_\mathrm{E}`$ is the Earth’s moment of inertia, $`m_\mathrm{E}`$ and $`m_\mathrm{M}`$ are the terrestrial and lunar masses, and $`G`$ is the gravitational constant.

As follows from *telescopic observations* of the Earth’s rotation with respect to the distant objects, collected over the last 300 years, $`\dot{T}_\mathrm{E}^{(\mathrm{tel})}=1.410^5`$ s/yr (Pertsev 2000). So, according to Eq. (3), $`\dot{R}^{(\mathrm{tel})}=2.53`$ cm/yr. On the other hand, immediate measurements of the Earth–Moon distance by *LLR technique* give an appreciably greater value $`\dot{R}^{(\mathrm{LLR})}=3.82`$ cm/yr (Dickey et al. 1994).

The most of attempts to explain the anomaly $`\mathrm{\Delta }\dot{R}=\dot{R}^{(\mathrm{LLR})}\dot{R}^{(\mathrm{tel})}=1.29`$ cm/yr were based on accounting for some additional geophysical effects (e.g. secular changes in the Earth’s moment of inertia $`I_\mathrm{E}`$). Unfortunately, they did not lead to a satisfactory quantitative agreement with observations.

Yet another promising interpretation of the above-mentioned anomaly $`\mathrm{\Delta }\dot{R}`$, from our point of view, is Hubble expansion in the local space environment, which should contribute to $`\dot{R}^{(\mathrm{LLR})}`$ but will not manifest itself in $`\dot{R}^{(\mathrm{tel})}`$. As follows from the standard relation $`\mathrm{\Delta }\dot{R}=H_0^{(\mathrm{loc})}R`$, the “local” Hubble constant should be

$$H_0^{(\mathrm{loc})}=\mathrm{\hspace{0.17em}33}\pm 5(\mathrm{km}/\mathrm{s})/\mathrm{Mpc}.$$
(5)
## 3 Discussion

At first sight, the Hubble constant given by Eq. (5) seems to be erroneous, since it equals only about one-half the commonly-accepted value at intergalactic scales. Nevertheless, it can be reasonably interpreted if the local Hubble expansion is formed only by some kind of an *unclumped* “dark matter” or “dark energy”, uniformly distributed in the Universe (such as $`\mathrm{\Lambda }`$-term, “quintessence”, inflaton-like scalar field, and so on); while the other kinds of matter experienced a gravitational instability, formed compact objects, and are no longer able to make their contributions to the rate of Hubble expansion at the local scales.

Besides, apart from any theoretical arguments, the recent observations (e.g. by Ekholm et al. 2001) revealed that a linear “quiescent” Hubble flow begins at least from the distances $``$1–2 Mpc, i.e. an order of magnitude less than was usually expected before. This fact was also interpreted by Chernin et al. (2001) as a manifestation of some unclumped dark matter.

So, the LLR technique, which was used so far only in solar-system studies, may also become a valuable tool for solving the cosmological problems, because it enables us either to measure a local rate of Hubble expansion or to impose an upper limit on this quantity (Dumin 2001a).

Of course, a cosmological nature of the anomalous increase in the Earth–Moon distance is not reliably established by now. It may be also caused merely by some geophysical artefacts. To distinguish between these two possibilities, it would be very interesting to seek a similar effect in some artificial satellite system, which does not suffer from geophysical uncertainties (Dumin 2001b). The suitable objects may be the space-based laser interferometers projected for searching the gravitational waves, such as LISA (e.g. Bender et al. 1998).

###### Acknowledgements.

I am very grateful to Yu.V. Baryshev, P.L. Bender, A.D. Chernin, E.V. Derishev, C. Hogan, S.M. Kopeikin, S.M. Molodensky, M. Tinto, and A.V. Toporensky for valuable discussions and critical comments, as well as to D.P. Kirilova and V.D. Kuznetsov for various help in the course of work.


text:
# New proposal for numerical simulations of 𝜃-vacuum like systems

(February 18, 2002)

## Abstract

We propose a new approach to perform numerical simulations of $`\theta `$-vacuum like systems, test it in two analytically solvable models, and apply it to $`\text{CP}^3`$. The main new ingredient in our approach is the method used to compute the probability distribution function of the topological charge at $`\theta =0`$. We do not get unphysical phase transitions (flattening behavior of the free energy density) and reproduce the exact analytical results for the order parameter in the whole $`\theta `$-range within a few percent.

Quantum Field Theories with a topological term in the action are a subject of interest in high energy particle physics and in solid state physics since a long time. In particle physics, these models describe particle interactions with a CP violating term. The extremely small experimental bound for the CP violating effects in QCD (strong CP problem) is still waiting for a convincing theoretical explanation strongcp . In solid state physics, chains of half-integer quantum spins with antiferromagnetic interactions are related to the two-dimensional $`O(3)`$ nonlinear sigma model with topological term at $`\theta =\pi `$. It has been argued that this model presents a second order phase transition at $`\theta =\pi `$, keeping its ground state CP symmetric (Haldane conjecture) haldane .

Non-perturbative studies of field theories with a $`\theta `$-vacuum term are enormously delayed because of the complex character of the euclidean action which forbids the application of all standard Monte-Carlo algorithms. Besides this, lattice QCD lacks from a simple consistent definition of topological charge.

The partition function $`𝒵_V(\theta )`$ of any $`\theta `$-vacuum like model in a finite space-time lattice volume $`V`$ can be written, up to a normalization constant, as the discrete Fourier transform of the probability distribution function (p.d.f.) of the topological charge at $`\theta =0`$:

$$𝒵_V(\theta )=\underset{n}{}p_V(n)\mathrm{e}^{\mathrm{i}\theta n},$$
(1)
where $`p_V(n)`$ is the probability of the topological sector $`n`$ and the sum runs over all integers $`n`$. In almost all practical cases the sum in (1) has a number of terms of order $`V`$ since the maximum value of the topological charge at finite volume is of this order.

Since efficient algorithms for numerical simulations of physical systems with complex actions are not yet available, the only a priori reliable numerical scheme to analyze the thermodynamics of $`\theta `$-vacuum like models goes through the determination of the p.d.f. of the topological charge, $`p_V(n)`$, and the evaluation of its Fourier transform (1). But this is a difficult task due to the following two technical reasons, which will be clarified later: i) any numerical determination of $`p_V(n)`$ suffers from statistical fluctuations PRD , and small errors in $`p_V(n)`$ can induce enormous relative errors in the determination of a quantity as $`𝒵_V(\theta )`$ which is an extremely small number of order $`e^V`$, ii) even if we were able to evaluate $`p_V(n)`$ with infinite precision, the sum in (1) contains terms that differ by many orders of magnitude, running from 1 to $`e^V`$ RANDOM .

In few specific cases one can overcome the sign problem STURN . However previous attempts by other groups to simulate $`\theta `$-vacuum like systems GENER ; PRD , were based on the numerical determination of the p.d.f. of the topological charge straightforwardly, by standard simulations, or by more sophisticated methods based on the use of multi-binning and re-weighting techniques. In all these attempts, artificial phase transitions at a $`\theta _c`$ decreasing with the lattice volume were observed for the two-dimensional U(1) gauge theory at strong coupling as well as $`\text{CP}^N`$ models. The origin of this artificial behavior, which follows from a flattening behavior of the free energy for $`\theta `$-values larger than a certain $`\theta _c`$, was analyzed in PRD ; JAPAN . Both groups agreed that the observed behavior was produced by the small statistical errors in the determination of the p.d.f. of the topological charge, the effect of which became more and more relevant as the lattice volume was increased. In Ref. JAPAN  it was also noticed that by smoothing the p.d.f. flattening disappears.

The purpose of this paper is to introduce a new numerical approach to simulate $`\theta `$-vacuum like models. This approach is based on a new method to compute the p.d.f. of the topological charge and the use of a multi-precision algorithm in order to compute the sum in (1) with a precision as high as desired.

For reasons which will become apparent in what follows, let us write the partition function (1) as a sum over the density of topological charge $`x_n=n/V`$ and set $`p_V(n)=\mathrm{exp}[Vf_V(x_n)]`$, where $`f_V(x)`$ is a smooth interpolation of $`1/V\mathrm{ln}p_V(x_n)`$:

$$𝒵_V(\theta )=\underset{x_n}{}e^{Vf_V(x_n)}\mathrm{e}^{\mathrm{i}\theta Vx_n}.$$
(2)
Equation (2) defines a $`2\pi `$ periodic function of $`\theta `$. Since CP is a symmetry of the action at $`\theta =0`$ and $`\theta =\pi `$, $`f_V(x)`$ will be an even function. We will assume that CP is realized in the vacuum at $`\theta =0`$ since otherwise the theory would be ill-defined at $`\theta 0`$ NOS . This implies that $`\mathrm{exp}[Vf_V(x_n)]`$ will approach a delta distribution centered at the origin in the infinite volume limit. In some exceptional cases, as QCD in the chiral limit, the function $`f_V(x)`$ is not defined since any topological sector with non vanishing charge has a vanishing probability. However, this is a trivial case in which the theory is independent of $`\theta `$.

Let us consider the partition function (2) in the complex $`\theta `$ plane, in particular on the imaginary axis $`\theta =ih`$, and let $`f(x)`$ be the infinite volume limit of $`f_V(x_n)`$. All the terms entering equation (2) are positive definite, and then in the infinite volume limit the free energy is given by the saddle point. Assuming that $`f(x)`$ has first derivative for any $`x`$ except at most in isolated points, we can write the saddle point equation:

$$f^{}(x)=h$$
(3)
which gives the external ”magnetic field” h as a function of the density of topological charge $`x`$.

Our proposal to compute the function $`f(x)`$ is based in the following three steps:

* To perform standard numerical simulations of our system at imaginary $`\theta =ih`$ and to measure the mean value of the density of topological charge as a function of $`h`$ with high accuracy (tipically a fraction of percent). This is feasible since the system we have to simulate has a real action. Then, Eq. (3) is used to get a numerical evaluation of $`f^{}(x)`$.
* To get $`f(x)`$ we have to integrate $`f^{}(x)`$. Between the possible ways to do this integral, we decided to fit $`f^{}(x)`$ by the ratio of two polynomials, whose order is chosen to obtain a high quality fit, and then to perform analytically the integral of the fitting function. In this way we get a very precise determination of $`f(x)`$, which allows us to compute the p.d.f. in a range varying several thousands orders of magnitude. This is the main advantage of our approach when compared with other methods based on a direct computation of $`p_V(n)`$.
* To use a multi-precision algorithm to compute the partition function (2) using as input the function $`f(x)`$ previously determined.

The function $`f(x)`$ obtained in step ii) suffers from statistical and systematic errors, the last coming from the fact that the saddle point equation (3) has finite volume corrections. An analysis of these errors for the models and sizes we have studied (see below), shows that systematic errors due to finite volume effects are smaller than the statistical ones in the whole relevant range of $`x`$. This is the reason why we will replace $`f_V(x_n)`$ in equation (2) by its asymptotic value $`f(x_n)`$ in what follows. This substitution has no effect in the infinite volume limit at imaginary $`\theta `$ and we are assuming that the same holds for real $`\theta `$.

Before presenting our results for the various testing models we have analyzed, let us briefly discuss how errors in the determination of $`f(x)`$ can propagate to $`𝒵_V(\theta )`$. This is an important point since, as pointed out in PRD ; JAPAN , the artificial phase transitions found in $`U(1)`$ and $`\text{CP}^N`$ were caused by the statistical errors in the determination of the p.d.f.

To this end, let $`f_V(x_n)`$ be the exact value at a given lattice volume $`V`$ of the function which parameterizes the p.d.f. and be $`\mathrm{\Delta }f_V(x_n)`$ a given deviation from the exact value. If we denote by $`𝒵_{}^{}{}_{V}{}^{}(\theta )`$ the partition function computed with $`f_V(x_n)+\mathrm{\Delta }f_V(x_n)`$, a simple calculation tells us that this partition function is related with the exact partition function $`𝒵_V(\theta )`$ as follows:

$$𝒵_{}^{}{}_{V}{}^{}(\theta )=𝒵_V(\theta )e^{V\mathrm{\Delta }f_V(x_n)}.$$
(4)
We can at this point analyze two extreme cases. First let us assume that $`\mathrm{\Delta }f_V(x_n)`$ vanishes everywhere except at a given $`x_m`$. Taking into account that the partition function $`𝒵_V(\theta )`$ which enters in the denominator of the expectation value $`e^{V\mathrm{\Delta }f_V(x_n)}`$ should behave as $`e^{Vg_V(\theta )}`$, where $`g_V(\theta )`$ is the free energy density ($`g_V(0)=0`$), we should get

$`e^{V\mathrm{\Delta }f_V(x_n)}=\mathrm{\hspace{0.33em}1}`$
$`+\mathrm{\hspace{0.25em}2}e^{V(f_V(\frac{m}{V})g_V(\theta ))}(e^{V\mathrm{\Delta }f_V(\frac{m}{V})}1)\mathrm{cos}(m\theta ).`$ (5)
Since $`g_V(\theta )`$ is an increasing function of $`\theta `$, it will be smaller than $`f_V(x_m)`$ near the origin, and therefore the free energy density computed with the modified partition function will differ from the exact one of a small quantity of order $`\frac{1}{V}`$ in this region. However at larger values of $`\theta `$ the function $`g_V(\theta )`$ can become larger than $`f_V(x_m)`$ and in such a case the correction to the free energy density will be finite or, even worst, the modified partition function can become negative. The other extreme case is that in which we assume that the error $`\mathrm{\Delta }f_V(x_n)`$ is constant. Under such assumption Eq. (4) implies that the error in the free energy density $`g_V(\theta )`$ will also be $`\mathrm{\Delta }f_V`$.

The previous discussion, as the results of JAPAN  suggests that correlated errors propagate in a less dramatic way than uncorrelated ones, and this would be a good scenario for numerical methods which, as in our approach, produce correlated errors in the determination of $`f_V(x_n)`$. We have checked this interesting issue in the two-dimensional U(1) model at strong coupling (see below) and verified that in fact correlated errors in $`f_V(x_n)`$ induce errors of the same order in the free energy density $`g_V(\theta )`$.

To test these ideas we have analyzed three models: the one-dimensional antiferromagnetic Ising model within an external imaginary magnetic field, the two-dimensional compact U(1) model with topological charge, and $`\text{CP}^3`$ in two dimensions. The coupling to the imaginary magnetic field in the Ising model can be written as $`\mathrm{i}\theta k_BT/2_iS_i`$, where $`k_B`$ is the Boltzmann constant and $`T`$ the temperature. For an even number of spins, the quantity $`1/2_iS_i`$ is an integer taking all values between $`N/2`$ and $`N/2`$, and therefore it can be seen as a quantized charge. Furthermore, the theory has a $`Z_2`$ symmetry at $`\theta =0`$ and $`\theta =\pi `$ which is the analogous of CP in field theory. We use the notation $`F=J/k_BT`$, where $`J`$ is the coupling constant between nearest neighbors. The transfer matrix technique allows to solve analytically the model. For antiferromagnetic couplings, $`F<0`$, the magnetization is an analytic function of $`\theta `$ between $`\pi `$ and $`\pi `$. At $`\theta =\pi `$ the system shows a first order phase transition with a non-vanishing magnetization. From a numerical point of view the determination of the free energy density and order parameter through equation (1) in this model has the same level of complexity of more sophisticated models. Furthermore, in contrast to two-dimensional $`U(1)`$ gauge theory, where the p.d.f. of the topological charge is nearly gaussian, the non-gaussian behavior of the p.d.f. of the mean magnetization in the antiferromagnetic Ising model makes this model a good laboratory to check the reliability of our approach. Fig. 1 shows our numerical results for the order parameter versus $`\theta `$ for a linear chain of 1000 spins and $`F=1/2`$. Statistical errors were estimated by doing 10 samples of the numerical results and applying a Jack-Knife analysis. The p.d.f. of the order parameter for such a system takes values in a range of around 2000 orders of magnitude. Notwithstanding that, we are able to reproduce the order parameter in the whole $`\theta `$ interval within a few percent.

The two-dimensional compact $`U(1)`$ gauge model with $`\theta `$-angle at strong coupling constitutes another interesting check because we can compare the goodness of our approach with the other existing simulations which showed artificial behavior with a fictitious phase transition moving to the origin when increasing the lattice volume. Fig. 2 displays our results for the topological charge density versus $`\theta `$ in a 80$`\times `$80 lattice at $`\beta =0`$ and $`\beta =0.6`$. We are able to reproduce the exact result within a few per thousand in the whole $`\theta `$ interval. The agreement between analytical and numerical results is actually impressive. Furthermore the flattening found in PRD  for the free energy density in relatively small lattices is absent in our simulations even in the 80$`\times `$80 lattice.

To test how different kind of errors in the determination of the function $`f(x)`$ which defines the p.d.f. of the density of topological charge can affect the determination of the free energy and order parameter, we have added to the measured $`f(x)`$ a random relative error of order $`10^3`$. Fig. 3 shows the order parameter obtained in this way. As can be seen a small but random error in $`f(x)`$ propagates to the order parameter in a very dramatic way and makes the calculation meaningless. Contrary to that if, in order to simulate a correlated relative error of order up to 50%, we replace the measured $`f(x)`$ by the (even) function $`f(x)\left(1+0.5\mathrm{sin}(x^2)\right)`$, the result for the order parameter is practically indistinguishable from the exact value for $`\theta <\pi /2`$, and the maximum deviation is about 25%, at $`\theta =\pi `$ (see Fig. 3). We conclude that random errors in $`f(x)`$ propagate in a very dramatic way but correlated errors do not, and this helps one to understand why our approach works so well.

The last model we have analyzed is $`\text{CP}^3`$ in two dimensional euclidean space. It is the standard wisdom that this model shares many qualitative features with QCD. Even if it has not been analytically solved we believe it is worthwhile to compare our results with previous existing numerical simulations. We studied the lattice formulation that makes use of an auxiliary U(1) field. Also in this model, the previous numerical simulations gave artificial phase transitions with a flattening behavior for the free energy density at a $`\theta _c`$ decreasing with the lattice volume GENER ; PRD .

Fig. 4 shows our results for the order parameter versus $`\theta `$ at $`\beta =0.6`$ on a $`100^2`$ lattice. We have chosen this particular $`\beta `$ value in order to compare directly with results reported in PRD . In our simulations we have no trace of the fictitious phase transition found in PRD . Furthermore the order parameter is clearly different from zero at $`\theta =\pi `$, hence CP is spontaneously broken at this $`\beta `$ value. An open question is how CP is realized in the continuum limit BURK . A more extensive analysis of this model, including a study of the behavior of the order parameter in the continuum limit, and an analysis of statistical and systematic errors involved in our approach, will be published elsewhere. What is interesting to notice here is that our results do not suffer from artificial phase transitions, caused by statistical errors in the determination of the p.d.f.

In the three models analyzed the finite size effects in $`f^{}(x)`$ cannot be appreciated since they are completely masked by the small statistical errors. For instance, finite size effects can be exactly computed in the Ising model: they are exponentially small with the lattice size. This is a general feature of non-critical systems. However, volume effects might be troublesome in the analysis of the continuum limit. Concerning systematic errors due to the choice of a particular fitting function for $`f^{}(x)`$, the difference between the numerical and exact results for the Ising and compact $`U(1)`$ models (beside the statistical errors) reported in Figs. 1 and 2 give us an idea on the order of magnitude of these errors. Of course systematic errors can depend on the model as well as on the parameters. In $`CP^3`$ at $`\beta =0.6`$ we did also a 5-parameters polynomial fit of the data for $`f^{}(x)`$ in the relevant $`x`$-interval. The discrepancy between the topological charge density obtained in the two cases is at most $`7\%`$ (see Fig. 4).

Similar ideas to those presented in this work have been proposed and promisingly applied to a matrix model of QCD at finite density in NISH .

This work has been partially supported by an INFN-CICyT collaboration and MCYT (Spain), grant FPA2000-1352. Victor Laliena has been supported by Ministerio de Ciencia y Tecnología (Spain) under the Ramón y Cajal program. Vicente Azcoiti and Victor Laliena would like to thank Professor Domingo Gonzalez for his invaluable contributions and very kind support.


text:
# A first principles study of sub-monolayer Ge on Si(001)

## 1 Introduction

Understanding the heteroepitaxial growth of Ge on Si(001) is vitally important for two reasons: first, it is a prototypical system for strained, Stranski-Krastanow growth; second, it has great potential for growing new semiconductor devices while remaining compatible with existing Group IV technology. We present here a first-principles investigation of the energetics of the $`(2\times n)`$ reconstruction observed at low Ge coverage, with the aim of explaining why the periodicity $`n`$ has the values observed experimentally; we also present an interpretation of the energetics in terms of the relaxed atomic geometries.

The clean Si (001) surface shows the well-known reconstruction due to the formation of rows of Si dimers. During the early stages of solid-source molecular beam epitaxy (SSMBE) of Ge on Si (001), the growth is remarkably similar to that of Si itself, with rows of dimers forming on the surface. As the coverage approaches one monolayer, however, the system shows the effects of strain due to the mismatch of the Si and Ge lattice parameters, and trenches of missing dimers appear (see Fig. 1). These trenches are oriented at right angles to the axis of the dimer rows, and there is a fairly regular spacing $`n`$ between trenches . This is the $`(2\times n)`$ reconstruction, with $`n`$ observed to be about 8, though the distribution of $`n`$ is fairly broad and somewhat dependent on growth conditions, with values up to 12 being reported. A similar reconstruction is observed during gas-source MBE (GSMBE) , though there are small differences due to the presence of hydrogen on the surface.

As further Ge is deposited (both in SSMBE and GSMBE), a series of further reconstructions are seen, which vary depending on growth conditions. The temperature at which growth occurs, the deposition rate of Ge and the presence of hydrogen all affect the growth . However, the next stage observed after the formation of $`(2\times n)`$ is generally the $`(m\times n)`$ reconstruction . This forms when a further layer of Ge grows on a $`(2\times n)`$ surface, and does not fill in the trenches in the surface. The increased strain in the new layer generally means that the new periodicity ($`m`$) is less than the old ($`n`$). On Si (001), there are two step types: parallel to the dimer rows (A-type, generally smooth) and perpendicular to the dimer rows (B-type, generally rough) . After one or two monolayers of Ge have been deposited, this roughness is seen to first equalise and then reverse . Finally, large scale features form along the elastically soft (100) and (010) directions: “hut” pits and clusters . Although the work reported here focuses only on the $`(2\times n)`$ reconstruction, it provides the foundation for investigating the other more complicated reconstructions. In particular, we will report elsewhere  on tight-binding investigations of the $`(m\times n)`$ reconstruction, which use the present first-principles calculations to validate a tight-binding parameterisation of the Si/Ge system for this kind of application.

The calculations to be presented are based on density-functional theory (DFT) in the generalized-gradient approximation (GGA), implemented with pseudopotentials and plane-wave basis sets. Starting from the perfect Ge monolayer as point of reference, we have calculated the formation energy and relaxed structure of missing-dimer trenches having spacings $`n`$ ranging from 4 to 12. By calculating the trench formation energy both with and without atomic relaxation, we shall show that the strong dependence of formation energy on $`n`$ arises almost entirely from relaxation effects. We shall also show that, in a sense that is appropriate to the usual experimental conditions, this formation energy is a minimum for a value of $`n`$ in the experimentally observed range.

There has been previous modelling of the Si/Ge $`(2\times n)`$ reconstruction, but only using empirical potentials, such as the Stillinger-Weber and modified Keating forms . These also studied atomic relaxation and strain effects, relating these to the periodicity of the reconstruction, and found results in broad agreement with experiment. Our calculations provide support for the physical mechanisms that emerged from these empirical studies. One problem addressed in the earlier work, but in our opinion not fully resolved, was that of comparing the energies of surface structures having different stoichiometries. This problem necessarily arises if one wishes to compare the energies of $`(2\times n)`$ reconstructions having different $`n`$ values, since removal or addition of Ge is needed to go from one to the other. The solution to this problem is equivalent to assigning an appropriate chemical potential to Ge, and different ways of doing this have been proposed. Since the predicted equilibrium value of $`n`$ depends on the choice of chemical potential, it is essential to identify the choice that corresponds to the real experimental conditions. We shall outline here what we believe to be the correct procedure. We note that during the experimental growth process there may be some intermixing of Ge and Si , though this can be suppressed by surfactants such as As or H. The likely effect is a reduction of surface strain, and hence an increase in the value of $`n`$; for the growth of the first monolayer, such effects are expected to be small, and we neglect them in the present work.

Technical details of our calculations are summarised in the next Section. We then present (Sec. 3.1) our results for energetics and relaxed geometry, first of the perfect Ge monolayer, and then of the missing-dimer trenches, including our analysis of the trench formation energy into electronic and relaxation contributions. Sec. 3.2 outlines our arguments about the correct Ge chemical potential to use in determining the equilibrium inter-trench spacing, and gives our numerical result for this, which is close to the experimental value $`n8`$. The paper ends with a summary of our conclusions.

## 2 Computational Methods

The fundamental ideas of DFT  have been extensively reviewed (see e.g. Refs. ), as have the pseudopotential and plane-wave techniques . The present calculations were performed using the VASP code , and employ the standard ultra-soft pseudopotentials  that form part of the code. The approximation we use for exchange-correlation energy is the generalised-gradient approximation (GGA) due to Perdew and Wang (PW91) . The choice of GGA rather than the local-density approximation (LDA) is deliberate. Since the energetics that interests us here depends quite sensitively on bonding and rebonding effects, and since the errors in bond energies are generally much larger with LDA than with GGA (LDA generally overbinds significantly), we regard the use of GGA as essential in this work.

DFT/pseudopotential/plane-wave calculations are most easily performed in periodic boundary conditions, and we therefore adopt the periodic slab geometry usually employed for surface-science work. The scientific issues we are addressing require the accurate treatment of quite small energy differences (typically on the order of 100 meV), and we have made efforts to ensure that the calculations are fully converged with respect to the thickness of the slabs and the width of the vacuum layer separating neighbouring slabs. For an eight layer $`(2\times 1)`$ reconstructed slab, we found that the total energy was converged to better than 1 meV for a vacuum layer of 5 Å (compared to 8 Å). For this vacuum width, we found that the change in surface energy in going from an eight layer slab to a twelve layer slab was less than 1 meV per dimer. (Both tests were conducted with a plane-wave cutoff energy $`E_{\mathrm{cut}}`$ of 225 eV and a $`4\times 4\times 1`$ Monkhorst-Pack k-point mesh.) Given the results of these tests, we chose to perform all the main Si/Ge calculations with eight-layer slabs, the top layer being Ge. Both top and bottom surfaces of the slab are reconstructed to form dimers, as further specified below. The width of the vacuum gap was taken to be 5 Å. Careful attention to basis-set completeness and Brillouin-zone sampling is also essential. In general, $`E_{\mathrm{cut}}`$ was chosen to be 225 eV, and $`k`$-point sampling was performed using the $`4\times 4\times 1`$ Monkhorst-Pack mesh, but detailed evidence will be presented below about the convergence of our results with respect to these parameters.

## 3 Results

### 3.1 Energetics and geometries

We begin by summarising the structure and energetics of the system having a perfect Ge monolayer in the $`p(2\times 2)`$ reconstruction (Fig. 2(a)). For the fully relaxed system, we find the following structural parameters (results from Refs.  in parentheses): Ge-Ge bond length 2.55 Å ($`2.382.44`$ Å); Ge-dimer tilt angle 19.2 ($`14.218.5^{}`$); up and down Ge-Si bond lengths 2.48 and 2.39 Å (2.42 and 2.34 Å). We note that the methods used to obtain the earlier results compared with here differ from ours in two significant ways: first, they used LDA rather than GGA; second, they used the $`(2\times 1)`$ rather than the $`p(2\times 2)`$ reconstruction. Both of these factors may have an appreciable effect on the structure of the monolayer; we note particularly that use of the $`(2\times 1)`$ reconstruction may well prevent the favourable relaxations along the dimer row allowed by $`p(2\times 2)`$. Given these differences, we regard the agreement with previous results as reasonable.

For comparison, we note the corresponding parameters for the $`p(2\times 2)`$ reconstruction on the clean Si (001) surface, obtained using the same GGA, and with the GGA lattice parameter for Si ($`a_0=5.45`$ Å, corresponding to a bulk bond length of 2.36 Å): dimer Si-Si bond length 2.36 Å; dimer-second layer up and down bond lengths 2.40 and 2.34 Å. The Ge-Ge dimer bond is elongated (actually beyond the bulk bond length, which we found to be 2.49 Å), with the extra freedom to relax allowing the substrate to take on more bulk-like lengths and angles.

We also report here the Ge monolayer formation energy per Ge dimer, denoted by $`E_\mathrm{m}`$, since it will be needed later. We define this to be the energy change per Ge dimer when we start with the perfect relaxed Si (001) $`2\times 2`$ surface and bring isolated Ge atoms from infinity to form the perfect relaxed $`2\times 2`$ Ge monolayer. We find the value $`E_\mathrm{m}=9.64`$ eV, and we have checked that this is converged within 0.01 eV with respect to slab thickness, plane-wave cut-off and $`k`$-point sampling.

We now turn to the energetics of formation of missing-dimer trenches. The results come from a series of calculations in which the periodically repeated cell contains a single missing-dimer trench, with spacings $`n=4`$, 6, 8, 10 and 12 between trenches. In all cases, the repeating cell is orthorhombic. Referring to Fig. 1 and denoting the lengths of the three cell edges by $`X`$ (along the dimer rows), $`Y`$ (along the missing-dimer trenches) and $`Z`$ (normal to the surface), these lengths are: $`X=(n/\sqrt{2})a_0`$, $`Y=\sqrt{2}a_0`$, and $`Z=15.88`$ Å, where $`a_0=5.44`$ Å. The value of $`Z`$ is appropriate to the eight-layer slab with vacuum width of 5.0 Å used in all the calculations (see Sec. 2).

We first report values for the fully relaxed missing-dimer formation energy $`E_\mathrm{f}(n)`$, defined as the energy per removed Ge-dimer needed to form an array of missing-dimer trenches with spacing $`n`$. In this process, we start from Si (001) with a non-defective monolayer of Ge, and the removed Ge atoms are taken to infinity, where they are isolated atoms in free space. In computational terms, we express $`E_\mathrm{f}(n)`$ as the difference $`E_{\mathrm{fin}}E_{\mathrm{init}}`$ of the fully relaxed energies per repeating cell of two systems, both having the same orthorhombic cell of dimensions $`(X,Y,Z)`$. The initial system (energy $`E_{\mathrm{init}}`$) is the slab with a single non-defective monolayer of Ge on one face. The final system (energy $`E_{\mathrm{fin}}`$) is formed from the initial system by removing one Ge-dimer from each repeating cell to create an array of infinite missing-dimer trenches with spacing $`n`$. We use exactly the same orthorhombic cells for the two systems, with the same plane-wave cut-off and $`k`$-point sampling, because this aids cancellation of errors.

We report in Table 1 our values for $`E_\mathrm{f}(n)`$ for the series of $`n`$ values. All the results were obtained for the eight-layer slab with the vacuum gap of 5 Å chosen for the reasons explained in Sec. 2. To show that our results are converged with respect to basis-set completeness and $`k`$-point sampling, we report values of $`E_\mathrm{f}(n)`$ for a moderate (150 eV) and a large (225 eV) plane-wave cut-off energy $`E_{\mathrm{cut}}`$ and for two $`k`$-point sets ($`4\times 4\times 1`$ and $`8\times 8\times 1`$). It is clear that $`E_\mathrm{f}(n)`$ is converged to within a few meV with respect to $`k`$-point sampling and plane-wave cut-off when we use a $`4\times 4\times 1`$ k-point mesh and a plane-wave cut-off of 225 eV. We note that $`E_\mathrm{f}(n)`$ is a monotonically decreasing function of $`n`$, which attains a plateau value for $`n>8`$. Its overall variation with $`n`$ is very substantial, since it is nearly 1 eV higher for $`n=4`$ than for $`n=12`$. This effective repulsion between missing-dimer trenches for spacings below $`n8`$ has been found before in calculations based on empirical models . Its consequence is that for a given overall density of trenches they will tend to become equally spaced.

Previous work  suggests that the effective repulsion between trenches stems from the elastic relaxation field surrounding each missing dimer. To test this, we have repeated the calculations of $`E_\mathrm{f}(n)`$, but without relaxation (we denote the unrelaxed value by $`E_\mathrm{f}^0(n)`$). As before, we start from the fully relaxed initial system, but when the Ge atoms are removed, all atoms are held fixed in their initial positions. The resulting $`E_\mathrm{f}^0(n)`$ values are reported as a function of $`n`$ in Fig. 3, and we see that their variation with $`n`$ is extremely small. This means that all the variation in the fully relaxed $`E_\mathrm{f}(n)`$ values comes from the relaxation of the final system with respect to the initial system. The relaxation energy $`E_\mathrm{f}E_\mathrm{f}^0`$ has a magnitude of nearly 2 eV for widely spaced missing dimers.

We show our calculated relaxed structure of the trench system in Fig. 2(b) for the spacing $`n=8`$. The key feature to notice is the large inward relaxation towards the trench, leading to Si–Si rebonding across the trench. Quantitatively, the relaxed Si–Si distance across the trench is 2.54 Å, to be compared with the bond length of 2.36 Å in the Si perfect crystal at ambient pressure. This is an elongation of $`8`$%, showing the strain that the system is under. It is interesting to compare this result with the single missing dimer in Si(001), which has an identical structure, apart from the Ge in the top layer. There, the bonds across the trench have a length of 2.56 Å ; the similarity suggests that the limiting factor on the relaxation in both cases is the Si-Si distance. The Ge dimers neighbouring the trench show an inward relaxation of 0.65 Å from their positions in the perfect monolayer. The up and down Si–Ge bond lengths neighbouring the trench are 2.58 and 2.45 Å, so that there is a significant lengthening compared with the perfect Ge monolayer. As expected, for smaller values of $`n`$, the relaxation around the trench is less pronounced. For example, in the $`n=4`$ case, the Si–Si separation across the trench is 2.63 Å, and the inward relaxation of the Ge dimers neighbouring the trench is 0.48 Å. The suppression of relaxation for small spacings is clearly responsible for the strong $`n`$-dependence of $`E_\mathrm{f}(n)`$, corresponding to the effective repulsion between trenches.

### 3.2 The experimental missing-dimer trench spacing

We now want to use our results for $`E_\mathrm{f}(n)`$ to clarify why the missing-dimer trenches adopt the spacings $`n`$ in the range $`812`$ observed in experiments. We explain first the statistical-mechanical basis for the arguments to be used. We assume that the surface is in thermal equilibrium, so that, for a given number of Ge atoms on the surface, the probability of finding any particular arrangement $`\gamma `$ of these atoms is proportional to $`\mathrm{exp}(E_\gamma /k_\mathrm{B}T)`$, where $`E_\gamma `$ is the energy of $`\gamma `$. (Strictly speaking, $`E_\gamma `$ should be a non-configurational free energy, but here we take it to be the equilibrium energy of $`\gamma `$.) Our thermal equilibrium assumption means that we are ignoring kinetic effects. We will discuss the validity of this assumption in Sec. 4.

Since our approach is to discuss the arrangements that will be seen in thermal equilibrium for a given number of Ge atoms on the surface, and we need only know how $`E_\gamma `$ varies as we go from one arrangement of these atoms to another, the energy zero chosen for $`E_\gamma `$ is irrelevant. However, it will be convenient to relate $`E_\gamma `$ values to the energy of one particular arrangement of Ge, which we call the ‘reference’ arrangement, whose energy is $`E_{\mathrm{ref}}`$. We choose this to be the arrangement in which all Ge atoms form perfect dimers, which are arranged to make a non-defective monolayer covering a certain area of the surface. The shape of the boundary of this area does not make any difference, but the following arguments become simpler if it is rectangular.

We now consider the energies of arrangements created from the reference arrangement by the formation of missing-dimer trenches. For definiteness, let there be $`L`$ dimer rows each containing $`P`$ dimers in the reference arrangement. We form $`Q`$ equally spaced missing-dimer trenches by removing $`QL`$ dimers, replacing them at the boundary of the monolayer, and allowing the whole system to relax. It is convenient to divide this process into two parts: (i) the prior fetching of $`2QL`$ Ge atoms from infinity and their deposition at the boundary of the reference system, in such a way that each dimer row is increased in length from $`P`$ dimers to $`P+Q`$ dimers, the number of dimer rows remaining the same; (ii) the subsequent formation of the missing-dimer trenches by the removal of $`QL`$ Ge dimers and their separation to infinity. We write the energy in process (i) as $`QLE_\mathrm{p}`$, where $`E_\mathrm{p}`$ is a constant energy that will be discussed further below. The energy change in process (ii) is $`QLE_\mathrm{f}(n)`$, where $`E_\mathrm{f}(n)`$ is the trench formation energy defined above, with the spacing $`n`$ given by $`n=(P+Q)/Q`$. The total energy change $`\mathrm{\Delta }E`$ is therefore:

$$\mathrm{\Delta }E=QL(E_\mathrm{f}(n)+E_\mathrm{p})=\frac{LP}{n1}(E_\mathrm{f}(n)+E_\mathrm{p})=LP\zeta (n),$$
(1)
where $`\zeta (n)=(E_\mathrm{f}(n)+E_\mathrm{p})/(n1)`$. The energetically most favorable value of $`n`$ is therefore obtained by minimizing $`\zeta (n)`$ with respect to $`n`$

In fact, the energy $`E_\mathrm{p}`$ of process (i) is the same as the monolayer formation energy $`E_\mathrm{m}`$ already discussed in Sec. 3.1. The reason is that the energy change on forming a certain amount of monolayer cannot depend on whether this is accomplished by placing Ge atoms on the clean Si surface, or by bringing Ge atoms to the boundary of a pre-existing piece of monolayer. (This assumes, of course, that we ignore edge effects, which is an appropriate approximation here.)

We now use our calculated $`E_\mathrm{m}=9.64`$ eV/dimer (see above, Sec. 3.1) to obtain numerical values for $`\zeta (n)`$, which we present in Fig. 4. We note that $`\zeta (n)`$ is negative for large $`n`$, so that it is energetically favourable to form widely spaced missing dimers. However, the repulsion between rows causes $`\zeta (n)`$ to increase at small $`n`$, and it has a minimum at $`n8`$. This optimum $`n`$ value corresponds well to the typical spacing observed experimentally.

## 4 Discussion and conclusions

A number of key points have emerged from our first-principles calculations on the $`(2\times n)`$ reconstruction of sub-monolayer Ge on Si (001). First, we have shown that the non-defective Ge monolayer is energetically unstable with respect to formation of widely spaced ($`n\mathrm{}`$) missing-dimer trenches. Second, there is a substantial effective repulsion between trenches, so that their formation becomes energetically unfavourable for small $`n`$. This means that in thermal equilibrium there is an optimal value of $`n`$, for which our calculations yield the estimate $`n8`$, in respectable agreement with the typical values observed experimentally. Third, we have shown that the effective repulsion is entirely due to elastic relaxation effects: the energy lowering due to relaxation is greatest when the strain fields of different trenches do not overlap, and decreases as $`n`$ decreases. Fourth, we have seen that for large $`n`$ there is significant rebonding between Si atoms in the trench, and this is probably essential in making trench formation energetically favourable.

The picture we have established is not entirely new. The work based on empirical models referred to in Sec. 1 came to essentially the same conclusion about the important role of elastic effects. The more fundamentally based first-principles calculations presented here therefore provide support for the earlier models. However, we have emphasised that a correct identification of the appropriate Ge chemical potential is crucial in understanding the energetic stabilisation that results from trench formation, as well as the equilibrium value of $`n`$. In earlier work, it was suggested that this chemical potential should be identified with the energy of bulk unstrained Ge (for large, fully relaxed islands), or else Ge biaxially strained to the Si lattice constant (for wide, coherent islands) . We have argued here that the formation energy of the perfect Ge monolayer provides the experimentally relevant point of reference in fixing this chemical potential.

Our reasoning is based on the assumption of full thermal equilibrium, and we have not attempted to account for kinetic effects. Since we are only trying to investigate the observed spacing $`n`$ of the reconstruction, this is reasonable. Although a variety of kinetic effects are observed during the formation of the $`(2\times n)`$ reconstruction (including the “displacive incorporation” growth model) , once the complete layer is formed it is stable, and kinetic effects are unlikely to play a role. We can also consider the limit of slow growth conditions, where thermal equilibrium will be a valid assumption . The observed range of $`n`$ is rather broad (roughly from 8 to 12), in part depending upon growth conditions; this is perfectly consistent with our results, which show only a weak increase in formation energy beyond the spacing $`n=8`$ at which it is a minimum. We also expect the temperature at which the growth occurs and the growth source to have an effect. We have also deliberately ignored intermixing between Ge and Si layers. There is MEIS evidence showing that for 1ML coverage (i.e. equivalent to the $`(2\times n)`$ surface we are modelling) at low temperatures (up to 500C) there is little intermixing , though recent measurements and calculations  indicate intermixing starting at about 500C. Nevertheless, the simplest model is one without intermixing, and this is where we have started. Intermixing will reduce the surface strain, thus increasing the value of $`n`$, as seen in experiment. Our value could thus be considered a lower limit, taken for thermal equilibrium and segregated Si and Ge layers.

JO and DRB thank EPSRC for funding through grants M01753 and M71640. DRB also thanks the Royal Society for funding through a University Research Fellowship. Calculations were performed on the Manchester CSAR service, through an allocation of time to the UKCP consortium (EPSRC grant M01753) and at the HiPerSPACE Centre at UCL (JREI grant JR98UCGI).