text:
# Quantum interference between multiple impurities in anisotropic superconductors

## Abstract

We perform a numerical study of the quantum interference between impurities in d-wave superconductors within a potential scattering formalism that easily applies to multiple impurities. The evolution of the low-energy local density of states for both magnetic and nonmagnetic point scatterers is studied as a function of the spatial configuration of the impurities. Further we discuss the influence of a subdominant bulk superconducting order parameters on the interference pattern from multiple impurities.

The past few years have proved the importance of experimental techniques which can directly test the wealth of information associated with modifications of the local density of states (LDOS) around impurities, grain boundaries and vortices in strongly correlated electron systems. In particular scanning tunneling microscopy (STM) measurements have provided detailed LDOS images around single nonmagnetic (Zn) and magnetic (Ni) impurities on the surface of the high temperature superconductor Bi2Sr2CaCuO4+δ (BSCCO).
For conventional superconductors Yu and Shiba showed that as a result of the interaction between a magnetic impurity and the spin density of the conduction electrons, a bound state located around the magnetic impurity is formed inside the gap in the strong-scattering (unitary) limit. For anisotropic superconductors a number of authors generalized the Yu-Shiba approach to study the LDOS around single impurities. It was found, for instance, that for a single nonmagnetic impurity there is only a virtual bound (or resonant) state due to the existence of the low-energy nodal quasi-particles. The one-impurity problem was recently reviewed by several authors.
Recently Hoffman et al. measured the energy dependence of the Fourier transformed LDOS images on the surface of optimally doped BSCCO below Tc. The dispersive features were explained from the point of elastic quasi-particle interference resulting from a single weak, nonmagnetic impurity. This gives credence that a scattering potential picture can yield valuable predictions in the superconducting state of these materials. Evidence for quantum interference between strong scatterers has been observed in the CuO chains of YBa2Cu3O6+x by Derro et al. Future experimental ability to control the position of the impurities on the surface of a superconductor and perform detailed STM measurements around multiple impurity configurations motivates further theoretical studies of the resulting quantum interference.
Previous calculations have studied the formation of bonding and antibonding states around two magnetic impurities in s-wave superconductors. For d-wave superconductors it was found that the interference effects between two nonmagnetic unitary scatterers depends crucially on the distance and orientation of the impurities. The orientational dependence arises from the anisotropic gap function, and provides an alternative method to identify the symmetry of the superconducting gap. Several authors have previously suggested similar ideas in the case of one impurity.
In this paper we study multiple impurity effects by exactly inverting the Gorkov-Dyson equation. In particular, we discuss the effect of quantum interference between: 1) nonmagnetic impurities in the strong scattering limit, 2) nonmagnetic impurities in the case of induced subdominant superconducting order parameters, and 3) magnetic and nonmagnetic impurities. All the calculations are performed within a quasi-particle scattering framework with classical impurities. The main purpose is to use quantum interference to obtain results that pose further tests on this approach and to illustrate the strong sensitivity of the LDOS on the positions of the impurities.
The Greens function $`\widehat{G}^{(0)}(𝐤,\omega )`$ for the unperturbed d-wave superconductor is given in Nambu space by

$$\widehat{G}^{(0)}(𝐤,\omega )=\left[i\omega _n\widehat{\tau }_0\xi (𝐤)\widehat{\tau }_3\mathrm{\Delta }(𝐤)\widehat{\tau }_+\mathrm{\Delta }^{}(𝐤)\widehat{\tau }_{}\right]^1,$$
(1)
where $`\widehat{\tau }_\nu `$ denotes the Pauli matrices in Nambu space, $`\widehat{\tau }_0`$ being the $`2\times 2`$ identity matrix and $`\widehat{\tau }_\pm =\widehat{\tau }_1\pm \widehat{\tau }_2`$. For a system with d$`_{x^2y^2}`$-wave pairing $`\mathrm{\Delta }(𝐤)=\frac{\mathrm{\Delta }_0}{2}\left(\mathrm{cos}(k_x)\mathrm{cos}(k_y)\right)`$. Below, $`\mathrm{\Delta }_0=25\text{meV}`$ and the lattice constant $`a_0`$ is set to unity. In this article we use a normal state quasi-particle energy $`\xi (𝐤)`$ relevant for BSCCO around optimal doping (14%)

$$\xi (𝐤)=2t\left(\mathrm{cos}(k_x)+\mathrm{cos}(k_y)\right)4t^{}\mathrm{cos}(k_x)\mathrm{cos}(k_y)\mu $$
(2)
with $`t=300\text{meV}`$, $`t^{}=0.4t`$ and $`\mu =1.18t`$. Here $`t`$ $`(t^{})`$ refers to the nearest (next-nearest) neighbor hopping integral and $`\mu `$ is the chemical potential.
We model the presence of scalar and magnetic impurities in the system by the following $`\delta `$-function interactions

$$\widehat{H}^{int}=\underset{i}{}\left[\left(V_i^S+V_i^M\right)\widehat{c}_i^{}\widehat{c}_i+\left(V_i^SV_i^M\right)\widehat{c}_i^{}\widehat{c}_i\right].$$
(3)
Here $`i`$ denotes the set of lattice sites containing magnetic and/or scalar impurities and $`V_i^M`$ ($`V_i^S`$) is the strength of the corresponding effective potential. We consider only the $`z`$-component of the magnetic impurity interaction and ignore spin-flip scattering.
For a single nonmagnetic impurity at the origin it is well-known that the scattering modifies the Greens function by

$$\delta G_{11}(𝐫,i\omega _n)=\frac{\left[G_{11}^{(0)}(𝐫,i\omega _n)\right]^2}{\frac{1}{V^S}G_{11}^{(0)}(0,i\omega _n)}\frac{\left[G_{12}^{(0)}(𝐫,i\omega _n)\right]^2}{\frac{1}{V^S}G_{11}^{(0)}(0,i\omega _n)}.$$
(4)
Here $`r`$ is the distance to the origin and $`G_{\alpha \beta }`$ refers to the $`\alpha \beta `$th entry of the $`2\times 2`$ Nambu subspace.
Naturally one can derive equivalent expressions for the LDOS modulations around several impurities. However, for a numerical study of the evolution of the LDOS for multiple impurities positioned at arbitrary lattice sites, we find it is easier to invert directly the real-space Gorkov-Dyson equation. The full Greens function $`\widehat{G}(𝐫,\omega )`$ is then obtained by solving the equation

$$\underset{¯}{\underset{¯}{\overset{^}{G}}}(\omega )=\underset{¯}{\underset{¯}{\overset{^}{G}}}^{(0)}(\omega )\left(\underset{¯}{\underset{¯}{\overset{^}{I}}}\underset{¯}{\underset{¯}{\overset{^}{H}}}^{int}\underset{¯}{\underset{¯}{\overset{^}{G}}}^{(0)}(\omega )\right)^1,$$
(5)
where the double lines indicate that the elements of this equation are matrices written in real- and Nambu space. The size of these matrices depends on the number of impurities and the dimension of the Nambu space. We have previously utilized this method to study the electronic structure around vortices that operate as pinning centers of surrounding stripes. We perform the 2D Fourier transform of the clean Greens function $`\widehat{G}^{(0)}(𝐤,\omega )`$ numerically by dividing the first Brillioun zone into a $`1400\times 1400`$ lattice and introducing a quasi-particle energy broadening of $`\delta =0.5\text{meV}`$ with $`\delta `$ defined from $`i\omega _n\omega +i\delta `$.
The differential tunneling conductance is proportional to the LDOS $`\rho (𝐫,\omega )`$ which in turn is determined from

$$\rho (𝐫,\omega )=\frac{1}{\pi }\text{Im}\left[G_{11}(𝐫,\omega )+G_{22}(𝐫,\omega )\right].$$
(6)
In the following we model the nonmagnetic unitary scatterers with a potential $`V^S=700\text{meV}`$ which gives rise to resonance energies around $`\pm 1.5\text{meV}`$ in agreement with experiment. (This is seen from the holelike resonance evident in the bottom LDOS scan in Fig. 5a. For a single impurity only one of the two resonances evident from Eqn. 4 has weight on the impurity site since the anomalous part of the Greens function, $`G_{12}(r,\omega )`$, vanishes at $`r=0`$ due to the symmetry of the d-wave gap.)
For interference between two nonmagnetic unitary impurities Morr et al. found strong variations in the LDOS as the distance between the impurities $`R`$ is varied along one of the crystal axes. The single-impurity spectrum was obtained when $`R`$ exceeds approximately $`10a_0`$. However, as expected for a $`d_{x^2y^2}`$-wave superconductor, this length scale is much larger along the nodal directions. This is seen from Fig. 1. Here the density of states is measured above one of the impurities fixed at the origin while the other is moved away along the nodal (a) or anti-nodal (b) direction. The single impurity LDOS is obtained for $`R`$ well above $`100a_0`$. Thus only for impurity concentrations below 0.1% does the LDOS correspond to the expected result from a single strong nonmagnetic impurity. For weaker scatterers the decay length will be considerably reduced.

For quantum interference between multiple fixed nonmagnetic unitary scatterers Fig. 2 shows the LDOS as the STM tip is scanned along a crystal axis on which the impurities are positioned.

In general the resonances are split by the proximity of other impurities and the number of resonances is directly proportional to the number of interfering impurities. However, locally the density of states may be strongly influenced by destructive interference. For instance, for the two impurities (Fig. 2a,c) sharp resonances exist only when $`R=2a_0`$ as is evident from Fig. 2c. When a third impurity is added at $`(2,0)`$ these resonances appear to broaden and shift to higher energies, Fig. 2d. Contrary to this, Fig. 2a,b show that the addition of a third impurity has only a minor effect when $`R=a_0`$.
The case of three nonmagnetic impurities is studied further in Fig. 3 which shows the evolution of the LDOS at $`(0,0)`$ as a function of the distance to a third impurity along the nodal (b,d) and anti-nodal (a,c) direction. The case without the third impurity corresponds to the topmost LDOS in Fig. 2a and 2c.

As in the case of two nonmagnetic impurities there are very strong variations in the final LDOS; the number of apparent resonances, their width and resonance energy depends crucially on the positions of the three impurities. The small modulations added by the third impurity seen in Fig. 3a,b agree with the destructive interference evident from the corresponding cases seen in Fig. 2a,b. Contrary to this, large modulations are again seen when increasing the distance between the two fixed impurities by a single lattice constant, Fig. 3c,d.
Recently Zhu et al. suggested a careful study of the two-impurity problem to extract information of the bulk Greens function of the clean system. In the following we show how the quantum interference between unitary scatterers is strongly affected by the induction of a small subdominant superconducting order parameter. Thus one may utilize the quantum interference between several impurities as an alternative method to detect a small subdominant order parameter.
For instance, tuning through a quantum phase transition from a d$`_{x^2y^2}`$ to a d$`_{x^2y^2}`$+idxy superconductor at a critical doping level, magnetic impurity concentration or magnetic field strength, a small dxy order could qualitatively alter the interference pattern. For $`\mathrm{\Delta }_{xy}(𝐤)=\mathrm{\Delta }_{xy}^0\mathrm{sin}(k_x)\mathrm{sin}(k_y)`$ with $`\mathrm{\Delta }_{xy}^0`$=5.0meV, we compare in Fig. 4 the LDOS at physically realizable impurity positions to the case with $`\mathrm{\Delta }_{xy}^0`$=0.

Also we show the difference between $`d+id`$ and $`d+is`$ pairing symmetry for these impurity configurations. For most spatial configurations the secondary pairing (id or is) leads to a sharpening of the resonances but at particular positions there is a qualitative difference as shown in Fig. 4. For instance, the induction of d+id pairing (Fig. 4b) can result in three apparent resonances contrary to the ground state with pure d$`_{x^2y^2}`$-wave pairing (Fig. 4a). Similarly, by comparing the LDOS at (1,1) (dashed lines) in Fig. 4d-f, it is evident that the interfering scatterers can provide a clear distinction between d+id and d+is pairing. Information of the induction of local order around the impurities can also be inferred from STM measurements of specific impurity configurations.
We turn now briefly to the study of the classical magnetic impurities in d$`_{x^2y^2}`$-wave superconductors. The interference of two magnetic scatterers in a s-wave superconductor was studied recently by Flatte et al. For comparison to the nonmagnetic case we use a magnetic potential strength $`|V^M|=700\text{meV}`$ which does, however, not model all magnetic impurities (e.g. Ni) on the surface of BSCCO. Future experiments will reveal whether the scattering potential formalism utilized here is appropriate or whether more correlated effects are required.
Fig. 5 shows the quantum interference between two unitary scatterers: (a,d) one magnetic and one nonmagnetic, (b,e) two parallel magnetic, and (c,f) two antiparallal magnetic. In all figures one impurity is fixed below the STM tip at the origin $`(0,0)`$ while the other is displaced along the: (a-c) horizontal crystal axis or (d-f) along the nodal direction. In Fig. 5a,d it is the nonmagnetic impurity that is fixed at the origin. Again the number of resonances, their position, amplitude and width depends in detail on the distance and nature of the two scatterers. Further, the spatial evolution of the LDOS is quite similar for case (a),(b), and (d),(e). These are, however, very different from the interference between two antiparallel magnetic impurities (c,f) which is dominated by strong destructive interference at small separations along the antinodal direction and a surprisingly fast recovery to the single impurity case along both the nodal and antinodal directions.

The results presented above remain qualitative since fits to a specific experiment in addition to details from the tunneling matrix elements could also include modified hopping integrals around the impurities, gap suppression and possibly both magnetic and nonmagnetic scattering potentials. We have checked that a gap suppression on the bonds around the impurity site does not produce any qualitative changes. However, on a phenomenological level the gap suppression could allow for a competing magnetic order parameter to develop around the impurity. Thus, the gap suppression might be important for explaining the formation of magnetic moments around nonmagnetic impurities as seen by NMR experiments. These issues are currently controversial but the vast amount of information inferred from the quantum interference between multiple impurities may help settle this, and more importantly settle the validity of the scattering potential scenario versus more correlated models.


format: text

text:
# Temporal distortion of annual modulation at low recoil energies

## I Introduction

It is widely believed, as suggested by a host of independent cosmological and astrophysical observations, that the most part of the matter in the Universe is not visible, revealing its existence only through gravitational effects. In particular, data on the rotational curves of galaxies indicate that the galactic visible parts are surrounded by approximately–spherical dark halos which extend up to several times the size of the luminous components.

The best candidates to provide dark matter in galaxies are Weakly Interacting Massive Particles (WIMP). Several WIMP direct–detection experiments are operating direct , with the goal of measuring the nuclear recoil energy (in the KeV range) expected to be deposited in solid, liquid or gaseous targets by the scattering of the non-relativistic dark halo WIMPs. Unfortunately, expected rates are small and the exponential decay of the WIMP recoil–spectrum resembles that of the background at low energies. However, a specific signature can be exploited in order to disentangle a WIMP signal from the background: the annual modulation of the rate modulation . This effect, expected to be of the order of a few per cent, is induced by the rotation of the Earth around the Sun. Due to its smallness, the annual modulation signature requires large–mass detectors with high statistics in order to overcome background fluctuations and be unambiguously detected. The annual modulation effect has been experimentally investigated by the DAMA Collaboration, which has indeed reported a positive evidence by using a 100 kg sodium iodide detector damalast .

One of the most important sources of uncertainty in the calculation of WIMP direct detection rates is the modeling of the velocity distribution function (DF) of the particles populating the dark halo. In the literature, a simple isothermal sphere model is usually adopted, i.e. a WIMP gas described by an isotropic Maxwellian with r.m.s velocity of the order of 300 Km s-1. This leads to a sinusoidal time–dependence of the expected signal with maximum (or minimum) around June 2nd, i.e. with the same (or opposite) phase as the relative velocity between the Earth and the halo rest frame.

However, the actual form of the WIMP velocity DF is unknown, and many different models, alternative to the isothermal sphere, are compatible with observations galaxy . The goal of the present Letter is to show that the main features of the annual modulation effect (the sinusoidal dependence with time, the occurrence of maxima and minima during the year and, under some circumstances, the one–year period) may be affected by anisotropies in the velocity DF. The most relevant effect is a distortion of the sinusoidal time–behaviour at low recoil energies. These energies, though below the current detector thresholds, might be reached in the future. The observation of the effects discussed in this Letter could provide informations on the phase space distribution of our galactic halo, especially on the degree of its anisotropy.

## II The annual modulation effect

Due to the rotation of the disk around the galactic center, the solar system moves through the WIMP halo, assumed to be at rest in the galactic rest frame. In the following, we will assume a right–handed system of orthogonal coordinates: the $`x`$ axis in the galactic plane, pointing radially outward; the $`y`$ axis in the galactic plane, pointing in the direction of the disk rotation; the $`z`$ axis directed upward, perpendicular to the galactic plane. Notice that our system differs from standard “galactic coordinates” by the different choice of the $`x`$ axis, which for us is directed outward.

The relative velocity between the WIMP halo and the detector is given by the Earth velocity $`\stackrel{}{v}^E`$, as seen in the galactic rest frame. It is the sum of three components: the galactic rotational velocity $`\stackrel{}{v}^G=(0,v_0,0)`$ Km s-1 (we will assume: $`v_0=220`$ Km s-1), the Sun proper motion $`\stackrel{}{v}^S=(9,12,7)`$ Km s-1langlsmith  and the Earth orbital motion $`\stackrel{}{u}^E(t)`$langlsmith :

$`v_x^E`$ $`=`$ $`v_x^G+v_x^S+u^E(\lambda )\mathrm{cos}\beta _x\mathrm{cos}[\omega (tt_x)],`$ (1)
$`v_y^E`$ $`=`$ $`v_y^G+v_y^S+u^E(\lambda )\mathrm{cos}\beta _y\mathrm{cos}[\omega (tt_y)],`$ (2)
$`v_z^E`$ $`=`$ $`v_z^G+v_z^S+u^E(\lambda )\mathrm{cos}\beta _z\mathrm{cos}[\omega (tt_z)],`$ (3)
where $`\lambda `$ is the ecliptic longitude, which is function of time. We can express $`\lambda `$ as langlsmith : $`\lambda =L+1^{}.915\mathrm{sin}g+0^{}.020\mathrm{sin}2g`$ where $`L=280^{}.460+0^{}.9856474t`$ and $`g=357^{}.528+0^{}.9856003t`$ and $`t`$ denotes the time expressed in days relative to UT noon on December 31. In Eqs.(13) $`u^E(\lambda )=u^E[1e\mathrm{sin}(\lambda \lambda _0)]`$ is the modulus of the Earth rotational velocity, which slightly changes with time due to the small ellipticity $`e`$ of the Earth orbit ($`u^E=29.79`$ Km s-1, $`e=0.016722`$ and $`\lambda _0=13^{}\pm 1^{}`$) langlsmith . In Eqs.(13) the $`\beta _i`$ denote the ecliptic latitudes and the $`t_i`$ are the phases of the three velocity components: $`\beta _x=174^{}.4697`$, $`\beta _y=59^{}.575`$, $`\beta _z=29^{}.812`$ and $`t_x=76.1`$ day, $`t_y=156.3`$ day, $`t_z=352.4`$ day. The angular velocity has a period of 1 year, and is given by: $`\omega =2\pi /(365\mathrm{days})`$. With the numbers given above, the modulus of the Earth velocity changes in time as: $`v_E|\stackrel{}{v}_E|233.5+14.4\mathrm{cos}[\omega (tt_0)]`$ (in Km s-1), where $`t_0152`$ days, i.e. June 2nd. Notice the slight offset between $`t_0`$ and $`t_y`$, due to the composition of the velocity components. As a good approximation, $`v_E`$ is usually taken as: $`v_E=(v_0+v_y^S)+u^E\mathrm{cos}\beta _y\mathrm{cos}[\omega (tt_0)]`$.

The direct detection differential rate $`dR/dE_R`$ is proportional to the integral:

$$(v_{\mathrm{min}})=_{wv_{\mathrm{min}}}𝑑\stackrel{}{w}\frac{f_{\mathrm{ES}}(\stackrel{}{w})}{w},$$
(4)
where $`f_{\mathrm{ES}}`$ and $`w`$ are the WIMP velocity DF and the WIMP velocity, respectively, in the Earth’s rest frame; $`v_{\mathrm{min}}`$ is the minimum value of $`w`$ for a given recoil energy $`E_R`$, WIMP mass $`m_W`$ and nuclear target mass $`m_N`$ and is given by: $`v_{\mathrm{min}}\sqrt{E_R/(2m_N)}(m_W+m_N)/m_W`$.

By indicating with $`f(\stackrel{}{v})`$ and $`\stackrel{}{v}`$ the WIMP velocity DF and the WIMP velocity in the galactic reference frame, the following transformations hold:

$`\stackrel{}{v}`$ $``$ $`\stackrel{}{w}=\stackrel{}{v}\stackrel{}{v}^E(t),`$ (5)
$`f(\stackrel{}{v})`$ $``$ $`f_{\mathrm{ES}}(\stackrel{}{w})=f(\stackrel{}{w}+\stackrel{}{v}^E(t)),`$ (6)
which imply that $`(v_{\mathrm{min}})`$, and so $`dR/dE_R`$, develops a time dependence induced by $`\stackrel{}{v}_E(t)`$.

## III The isotropic isothermal sphere

In the case of the isothermal sphere model the DF is given by a truncated isotropic Maxwellian which depends only on $`v|\stackrel{}{v}|`$:

$$f(v)=N\mathrm{exp}\left(\frac{v^2}{2\sigma ^2}\right)|_{vv_{\mathrm{esc}}},$$
(7)
where $`\sigma `$ is the WIMP r.m.s. velocity, given by: $`\sigma ^2=3v_0^2/2`$. It is clear that, through the change of reference frame of Eq. (6), $`(v_{\mathrm{min}})`$ depends on time only through $`v_E`$. Since the relative change of $`v_E`$ during the year is of the order of a few per cent, we can approximate $`(v_{\mathrm{min}})`$ with its first–order expansion in the small parameter $`ϵ\delta v_E/v_E`$, around its mean value $`_0(v_{\mathrm{min}})`$: $`(v_{\mathrm{min}})_0(v_{\mathrm{min}})\mathrm{cos}[\omega (tt_0)]`$. The well known result is then obtained that the WIMP rate has a sinusoidal time dependence with the same phase ($`t_0`$ June 2nd) as $`v_E`$, for all values of $`v_{\mathrm{min}}`$. This is shown in Fig.1, where $`(v_{\mathrm{min}})`$ is plotted as a function of time for various values of $`v_{\mathrm{min}}`$.

## IV Anisotropic models

The simplest generalization of the isothermal sphere model is given by a triaxial system described by a multivariate gaussian:

$$f(\stackrel{}{v})=N\mathrm{exp}\left(\frac{v_x^2}{2\sigma _x^2}\frac{v_y^2}{2\sigma _y^2}\frac{v_z^2}{2\sigma _z^2}\right)|_{vv_{\mathrm{esc}}},$$
(8)
where $`N`$ is the normalization constant. For $`v_{\mathrm{esc}}\mathrm{}`$, then $`N=[(2\pi \sigma _x^2)(2\pi \sigma _y^2)(2\pi \sigma _z^2)]^{1/2}`$. The usual isothermal sphere is the spherical limit of Eq.(8), obtained with: $`\sigma _x=\sigma _y=\sigma _z\sigma `$. In order to discuss the effect of anisotropy at fixed WIMP mean kinetic energy, we will fix $`\sigma ^2=\sigma _x^2+\sigma _y^2+\sigma _z^2`$ as in the isothermal case ($`\sigma ^2=3v_0^2/2`$) and discuss our results in terms of the two independent parameters: $`\lambda _{12}\sigma _x/\sigma _y`$ and $`\lambda _{32}\sigma _z/\sigma _y`$.

At variance with the isothermal sphere, now $`(v_{\mathrm{min}})`$ depends in general on all the three components of $`\stackrel{}{v}_E`$, and not simply on $`v_E`$. We can write:

$$f_{\mathrm{ES}}(\stackrel{}{w})=N\mathrm{exp}\left[(\stackrel{}{r}+\stackrel{}{\eta })^2\right],$$
(9)
where we have defined the reduced (dimensionless) variables: $`r_i=w_i/(\sqrt{2}\sigma _i)`$ and $`\eta _i=v_i^E/(\sqrt{2}\sigma _i)=\eta _i^0+\delta \eta _i\mathrm{cos}[\omega (tt_i)]`$ ($`i=x,y,z`$). In the isotropic case (isothermal sphere) one has $`\eta _{x,z}^01`$, $`\eta _y^01`$ and $`\delta \eta _i1`$. The presence of the order–one parameter $`\eta _y^0`$ and of the small oscillation amplitudes $`\delta \eta _i`$ allows the Taylor–expansion of $`(v_{\mathrm{min}})`$ in terms of $`ϵ_i\delta \eta _i/\eta _y^0`$ parameters. A straightforward calculation shows that the conclusions of the previous Section are recovered.

On the contrary, allowing anisotropies such that $`\lambda _{12}<1`$ and/or $`\lambda _{32}<1`$, the values of the parameters $`\eta _{x,z}`$ and $`\delta \eta _{x,z}`$ are enhanced, and the time–dependence of $`v_x^E`$ and $`v_z^E`$ in $`(v_{\mathrm{min}})`$ may become important. An example of this situation is shown in Fig. 2, where, the time evolution of $`(v_{\mathrm{min}})`$ is plotted for $`\lambda _{12}=0.2`$ and $`\lambda _{12}=0.8`$ (i.e.: $`(\sigma _x,\sigma _y,\sigma _z)=(42,208,166)`$ Km s-1.) This choice of $`\lambda _{12}`$, $`\lambda _{32}`$, which refers to a tangential anisotropy, corresponds to triaxial models discussed, for instance, in Refs. carollo ; galaxy . A distortion of the curves of Fig. 2, as compared to the familiar sinusoidal time–dependence, appears: this effect may be explained by the fact that now $`\delta \eta _x1`$, a Taylor–expansion of the type used in the isothermal sphere case breaks down and a full numerical calculation of the integral of Eq.(4) is required. The final result is not sinusoidal. This peculiar behaviour is more pronounced at low values of $`v_{\mathrm{min}}`$ (i.e. low recoil energies) namely for $`v_{\mathrm{min}}<`$ 80 Km s-1, for which the distortion is strong and the maxima (in absolute value) of the rate are shifted as compared to the standard case krauss . For larger values of $`v_{\mathrm{min}}`$ the distortion is less pronounced, and it dies away when $`v_{\mathrm{min}}>230`$ Km s-1. This may be explained by the fact that, as $`v_{\mathrm{min}}`$ grows, the integral of Eq.(4) becomes less sensitive to the parameter $`\eta _x`$ since it gets increasingly dominated by WIMPs with velocities along the $`y`$ axis, which is the one along which the boost due to the galactic rotation is directed. We notice that for values of $`v_{\mathrm{min}}`$ around ($`100200`$) Km s-1 a distortion is present, but the amplitude of the modulation is suppressed and therefore difficult to detect.

As a second example, in Fig.3 we plot $`(v_{\mathrm{min}})`$ as a function of time for the case $`\lambda _{12}=10`$, $`\lambda _{32}=3`$ (i.e.: $`(\sigma _x,\sigma _y,\sigma _z)=(257,26,77)`$ Km s-1. This situation is representative of a radial anisotropy. In this case, we have further enhanced the contribution of $`\delta \eta _y`$ over $`\delta \eta _{x,z}`$, so that one should expect to draw the same conclusions as in the isothermal sphere case, with the usual sinusoidal time–dependence of the rate and a phase close to $`t_yt_0`$. This is indeed the case, except for a narrow interval of values of $`v_{\mathrm{min}}`$ around $`v_{\mathrm{min}}`$ 210 Km s-1. In this range, $`(v_{\mathrm{min}})`$ develops two maxima because, for that particular choice of $`v_{\mathrm{min}}`$, there is an exact cancellation in the first term of the expansion in $`ϵ_y=\delta \eta _y/\eta _y^0`$, so that the term $`𝒪(ϵ_y^2)`$ proportional to $`\mathrm{cos}^2[\omega (tt_y)]`$ sets in. This particular cancellation of the first term in the Taylor expansion of $`(v_{\mathrm{min}})`$ happens also for the isothermal–sphere model, but in that case the size of the quadratic term is strongly suppressed because $`ϵ_y`$ is much smaller. In a WIMP direct detection experiment this effect would show up in a very peculiar way: a halving of the modulation period of the rate in a narrow range of recoil energies. Note, however, that in order to have some realistic chance to detect this effect, it should show up in one of the experimental energy bins just above threshold, where the highest signal/background ratio is usually attained.

In order to establish a link between our discussion and WIMP direct detection experiments direct , in Fig. 4 we plot $`v_{\mathrm{min}}`$ as a function of the quenched nuclear recoil energy $`E_{ee}qE_R`$ ($`q`$ is the quenching factor) for the target nuclei: Na, I, Ge, Te and for different WIMP masses. The vertical dashed lines show current energy thresholds achieved by each type of detector. Fig. 4 shows that values of $`v_{\mathrm{min}}<`$ 80 Km s-1, i.e. sufficiently low to observe a sizeable distortion effect as the one discussed for tangential anisotropy, correspond to WIMP recoil energies below the threshold of present direct detection experiments, and that the effect would be more easily detected at higher WIMP masses. However, a foreseeable lowering of the threshold, down to 0.5–1 KeV, could be enough to observe the distortion. On the other hand, for radial anisotropy, we can conclude that the recoil energy corresponding to a halving of the modulation period can actually coincide to the experimental thresholds within the reach of present–day detectors for 20 GeV $`<m_W<`$ 100 GeV, depending on the particular target nucleus. In this respect, we note that the properties of the annual modulation effect observed by the DAMA/NaI experiment damalast  (a one–year–period sinusoidal behaviour in the 2–6 KeV energy bins damalast ; galaxy ) implies that the DAMA/NaI experiment is already able to set constraint on strong radial anisotropies.

## V Conclusions

In the present Letter we have shown that the main features of the annual modulation of the signal of WIMP direct searches, i.e. the sinusoidal dependence of the rate with time, the position of its maxima and minima during the year and even the period, may be affected by relaxing the isothermal sphere hypothesis in the description of the WIMP velocity phase space. We have considered a multivariate gaussian and found that different situations may occur, depending on the pattern of anisotropy: tangential anisotropies induce a departure at low energies from the usual sinusoidal time–dependence, along with a shift in the position of the maximum of the signal during the year, while radial anisotropies may produce a halving of the modulation period in a particular energy bin. The former effect turns out to be relevant at low recoil energies, actually below the threshold of present–day experiments, while the latter should be already within the reach of current detectors. In particular, the properties of the annual modulation effect observed by the DAMA/NaI experiment damalast  may already indicate that strong radial anisotropies are excluded.

The effects discussed in this Letter may be used in the future to provide a direct way of measuring (or setting limits on) the degree of anisotropy of the galactic halo.


text:
# Auger Effect in the High-Resolution Ce 3d-edge Resonant Photoemission

## I INTRODUCTION

Investigations on the electronic structures of strongly correlated materials such as $`d`$ and $`f`$ electron systems have been very active lately. One of the perhaps oldest problems in this category is the origin of the $`\alpha \gamma `$ phase transition of Ce metal and its compounds, which has been studied for more than 50 years but still remains controversial. Photoemission spectroscopy, among various experimental techniques, is a very powerful tool that can directly probe the electronic structures of solids. Indeed in the case of Ce problem as well, the study of 4$`f`$ spectral weights of various Ce compounds using resonant photoelectron spectroscopy (RPES) technique at the Ce 4$`d`$-edge contributed decisively to the understanding of this famous phase transition. This technique is necessary to enhance the 4$`f`$ electron emissions and separate them from contributions of other conduction electrons, and utilizes the following process

$`4d^{10}4f^1+\mathrm{}\omega 4d^94f^24d^{10}4f^0+photoelectron.`$
The conventional wisdom that came out of these studies, at least until early 1990’s, was as follows; The major factor distinguishing between $`\alpha `$ and $`\gamma `$ phases of Ce and its compounds is the strength of the hybridization between Ce 4$`f`$ level and the conduction bands, and the disappearance of the local magnetic moment in the $`\alpha `$ phase results from the Kondo phenomenon of the singlet formation due to the hybridization of a localized 4$`f`$ electron with conduction electrons, and not from the promotion of the 4$`f`$ electron to the empty state above the Fermi level.

After the high resolution RPES was developed in the early 1990’s, however, some groups questioned this conventional wisdom and proposed a new interpretation of the 4$`f`$ spectral weights of Ce compounds. The main point of those groups was that the temperature or the material dependence of the 4$`f`$ spectral weights of Ce compounds near the Fermi level, which had been interpreted as the tail of the Kondo resonance, did not follow the Kondo scenario which stipulates that the Kondo temperature is the universal controlling parameter that determines the strength of the Kondo resonance near the Fermi level in different materials or at different measuring temperatures. This controversy generated many more careful studies on the 4$`f`$ spectral weights of various Ce compounds by several research groups,  some utilizing RPES technique and others using ordinary ultraviolet photoemission (UPS) technique with HeI (h$`\nu `$ = 21.2 eV) and HeII (h$`\nu `$ = 40.8 eV) sources.

Most of these later studies also supported Kondo resonance scenario, but one important issue which has not been settled so far is the surface effect and how to separate the surface from the bulk contributions. Since the electron kinetic energy at the Ce $`4d`$-edge RPES is about 120eV, and the electron mean free path at this kinetic energy is only about 5Å, the surface contribution in the 4$`f`$ spectral weights obtained with $`4d`$-edge RPES technique is expected to be very much significant, even larger than the bulk in some cases. The situation with the ordinary UPS data with HeI (h$`\nu `$ = 21.2 eV) and HeII (h$`\nu `$ = 40.8 eV) sources is not much different. Hence to compare quantitatively the experimental photoemission spectra with those expected from the Kondo temperature of bulk materials, it is imperative to be able to extract the bulk contributions from the experimental data. Many approaches have been tried with the Ce $`4d`$-edge RPES or ordinary UPS data, but the results were not conclusive enough since each method has its own drawbacks and uncertainties.

Recently a powerful new method to extract reliable bulk Ce 4$`f`$ spectral weights has become possible thanks to the advent of the new beamline at SPring-8, Japan. It has been well-known that we can reduce the surface contribution in photoemission spectra by using high energy photon source. For example, if around 880eV photon source is used in RPES, the electron mean free path of the photoelectrons from Ce 4$`f`$ levels will be about 20 Å. Since those 880eV photon energies can excite a Ce $`3d`$ core-electron to an empty 4$`f`$ level (Ce 3$`d`$ absorption edge), the 4$`f`$ spectral weight is also enhanced due to the resonance phenomenon as in the case of Ce 4$`d`$-edge RPES. Hence we can obtain much more bulk-sensitive Ce 4$`f`$ spectral weights by using those high photon energies in this so-called Ce 3$`d`$-edge RPES technique. In fact, this technique had been already tried more than 10 years ago, but the energy resolution in the Ce 3$`d`$-edge RPES at that time was not good enough to give conclusive results. Recently, Sekiyama et al. succeeded to obtain the total energy resolution including energy analyser at the Ce $`3d`$ edge photon source less than 100meV with high photon flux in the undulator beamline of SPring-8 and study high-resolution bulk-sensitive Ce spectral weights in CeRu2. This opened up a new possibility to extract reliable bulk-sensitive spectral weights in many interesting Ce compounds, and several papers have already appeared utilizing this Ce 3$`d`$-edge RPES technique.

However, in the course of this study, we realized that the experimental photoemission data depends sensitively on the exact incident photon energy near the Ce 3$`d`$-edge, and the difference of spectral lineshape is significant enough to affect data interpretation. Also the new Ce 3$`d`$-edge RPES data were found to be not completely consistent with the old data taken in early 1990’s, even after the difference of the energy resolutions is taken into account. To derive reliable conclusions from the Ce 3$`d`$-edge RPES data, therefore, it is necessary to understand the origin of these differences and changes. This paper is an attempt to understand the origin of these phenomena by systematically studying Ce 3$`d`$-edge RPES of several Ce compounds. We will show in this paper data on CeFe2, CeNi2, and CeSi2 only, although we studied many more Ce compounds which gave fully consistent results with the conclusion of this paper. The three Ce compounds to be discussed here have various Kondo temperature, the universal parameter which controls many physical properties of Ce compounds in the Kondo resonance scenario. The Kondo temperatures of intermetallic compounds CeFe2 and CeNi2 are larger than about 500 K, and that of CeSi2 compound is about 45 K. This paper will mainly focus on the general phenomena of the Ce 3$`d`$-edge RPES, and the electronic structures of individual compounds will be discussed in separate publications.

## II EXPERIMENTAL

Polycrystalline samples of CeFe2, CeNi2, and CeSi2 were made by the arc melting method under Ar gas environment with constituent elements Ce, Fe, Ni, and Si whose purity was better than 99.9 %. The homogeneity of samples was checked after annealing with x-ray diffraction.

Resonant photoemission spectra of the valence bands at various photon energies near the Ce $`3d`$-edge (h$`\nu `$ 880 eV) along with the Ce $`3d`$-edge x-ray absorption spectrum (XAS) were obtained at the beamline BL25SU in the SPring-8 of Japan. The base pressure of an analysis chamber was about 3 $`\times `$ 10-10 mbar or better. The samples were scraped with a diamond file to get clean surfaces, and every time the cleanness of sample surface was checked with the O $`1s`$ core-level spectrum. If the oxygen contamination peak was found noticeable in O $`1s`$ core-level spectrum, the samples were scraped again with the diamond file. All spectra were measured at the temperature of 20 K, and the scraping was also done at 20 K. The XAS measurements were done by the total electron yield mode, and the energy resolution was about 80 meV. Scienta 200 analyzer was used for the electron energy analysis for PES measurements, and the total resolution of the PES spectra was better than 100 meV.

## III DATA AND DISCUSSIONS

Fig. 1 shows the Ce MV XAS data representing $`3d_{5/2}4f`$ transition for CeFe2, CeNi2, and CeSi2 compounds. The $`3d_{5/2}`$ XAS of Fig. 1 consists of the main peak structure near $`h\nu `$ = 882 eV, and the satellite peak near $`h\nu `$ = 887 eV. The main peak structure comes from $`\mathrm{"}3d^{10}4f^1\mathrm{"}\mathrm{"}3d^94f^2\mathrm{"}`$ transition, and the satellite peak $`\mathrm{"}3d^{10}4f^0\mathrm{"}\mathrm{"}3d^9f^1\mathrm{"}`$ transition. The mark of ” ” means approximate interpretation, since the complete eigenstate of each peak is composed of many body configurations. The lineshape of the main peak near 882 eV is primarily determined by the multiplet structures of the $`\underset{¯}{3d}4f^2`$ electronic configuration, where the underline represents a hole. At the bottom of the figure is shown the theoretically expected XAS line spectrum, which is calculated by considering all the multiplets of the $`\underset{¯}{3d_{5/2}}4f^2`$ electronic configuration along with the appropriate transition probability from the ground state. We can see that the theoretical calculation compares favorably with the experimental XAS data, especially for the case of CeSi2, although we can also notice that the details of the main peak lineshape and the satellite peak intensity near 887eV change somewhat depending on the compounds. We will deal with this phenomenon and its implication on the electronic structures systematically in a separate publication. In this figure, the symbols A to G designate the incident photon energies used to obtain the $`4f`$ spectral weights with the Ce $`3d`$-edge RPES to be shown in subsequent figures.

Fig. 2 represents RPES spectra near the Fermi level of CeFe2 obtained with various incident photon energies corresponding to the symbols from A to G in Fig. 1. All through these photon energies the Ce 4$`f`$ spectral weights are resonantly enhanced by the following process.

$`3d^{10}4f^1+\mathrm{}\omega 3d^94f^23d^{10}4f^0+photoelectron.`$
Since the contribution from the Ce 4$`f`$ emissions is at least ten times larger than those from other conduction electrons at these resonance energies, we can safely neglect the other conduction electron contributions and consider the spectra of fig. 2 as the Ce 4$`f`$ spectral weights near the Fermi level, as in the case of Ce 4$`d`$-edge RPES. However, the bulk contribution of Ce 4$`f`$ spectral weights of Fig. 2 is much larger than that of $`4d`$ edge RPES, since the electron mean free path at $`3d4f`$ transition energies is much longer.

The spectra in Fig. 2 consist of a peak near the Fermi level and a shoulder at about 2eV. The shoulder near 2eV comes from the $`f^1f^0`$ transition, whereas the peak near the Fermi level comes from $`f^1f^0\underset{¯}{c}`$ transition(c means the hole in the conduction band) which corresponds to the tail of the Kondo resonance peak. The detailed analysis on those peaks and their relation to the electronic structures of CeFe2 have been dealt with in other papers. Here we focus on the change of the 4$`f`$ spectral lineshapes depending on the incident photon energy. We clearly see that as the photon energy is increased from A (881.76eV) to F (882.48eV), which is only 0.72eV apart, the width of the peak near the Fermi level is significantly broadened. To see these changes of spectral weights more quantitatively, we subtracted the spectrum at A from data taken at higher photon energies after normalizing each spectrum at its maximum intensity. The resulting difference spectra at various photon energies are shown in Fig.3.

We find two noticeable features in every difference spectrum of Fig.3 — a broad peak below the Fermi level EF and a sharp dip at EF. The maximum of the broad peak moves to the high binding energy side as the incident photon energy is increased. In fact, we note that the binding energy of the maximum position is exactly the same as the difference of two incident energies, which is indicated by the arrows. This tells us the kinetic energy of the broad peak remains the same regardless of the incident photon energy, a telltale signature of Auger emission. Hence, we can conclude that this broad peak results from the Auger transition, and that the apparent change of the spectral lineshapes depending on the incident photon energy is due to these Auger emissions overlapping the resonant 4$`f`$ photoemissions. This Auger process must be the Coster-Kronig type between different multiplets of the $`3d^94f^2`$ configuration in the intermediate state of RPES.

To obtain genuine Ce 4$`f`$ spectral weights from RPES data, therefore, it is necessary to avoid the contribution from these overlapping Auger emissions. The best way to achieve this is probably to take the spectrum at the incident photon energy corresponding to the lowest multiplet of the intermediate state. Indeed, we can see in the figure that the difference spectrum between B and A is nearly flat, implying that the spectral shape is nearly identical for data taken near the lowest multiplet. We have also confirmed that the spectra taken at incident photon energies slightly below A remain identical in lineshapes as that of A. This fact gives us an important lesson in extracting bulk-sensitive Ce 4$`f`$ spectral weights using the 3$`d`$-edge RPES technique — that is, it is dangerous to use the incident photon energy which simply gives maximum resonance effect (such as peak E in Fig.1), since at that energy the spectra might be contaminated with overlapping Auger emissions. It is instead best to choose the incident photon energy as low as possible which still gives appreciable resonance effect. Some of the previous Ce 3$`d`$-edge RPES experiments may have been performed at the photon energy of maximum resonance, in which case the data have to be interpreted with caution. This interpretation on the origin of the lineshape change with the incident photon energy can also explain the dip structure near EF shown in Fig.3. Because of the Auger transition between multiplets, the intensity of the Kondo resonance peak near EF will be reduced while some intensities move to higher binding energy side. Therefore if each spectrum is normalized to its maximum height before subtraction, as we did in our case, a dip will show up near the Fermi level with negative intensity in the difference spectra.

This phenomenon of overlapping Auger emissions in the Ce 3$`d`$-edge RPES spectra is not limited to the particular case of CeFe2. In Fig. 4 and Fig. 5 are shown the Ce 3$`d`$-edge RPES data for CeNi2 and CeSi2, respectively, taken at various incident photon energies. Again, we see that they both show changes of spectral shapes depending on the photon energy, similar to the case of CeFe2. The two spectra of CeNi2 in Fig.4 are obtained with the incident photon energy of 882.34eV (E) and 881.54eV (A), which are very close to each other and roughly corresponds to two maximum points of $`3d_{5/2}`$ XAS spectrum of CeNi2 in Fig. 1. But these two 4$`f`$ spectral weights have appreciably different shapes in that the spectrum at photon energy E shows a fairly strong peak around 0.8eV binding energy, which becomes very weak or absent in the spectrum A. To understand the origin of this change of spectral shape, we again took the difference spectrum of the two and showed it at the bottom of Fig. 4. We can see in this difference spectrum a broad peak with its maximum around 0.8eV binding energy, which is exactly the difference of the photon energy and therefore corresponds to the constant kinetic energy feature. This tells us that in this case of CeNi2 as well the origin of this phenomenon is the same as the case of CeFe2 discussed earlier, and the overlapping Auger emissions between different multiplets of $`3d^94f^2`$ intermediate state contribute to the lineshape change with the incident photon energy. The CeSi2 spectra shown in Fig. 5 gives the same story, although the phenomenon does not look as pronounced as in the case of CeNi2.

## IV CONCLUSION

In the newly-developed Ce 3$`d`$-edge resonance photoemission technique for the study of bulk electronic structures of Ce compounds, the spectral shapes were found to depend sensitively on the incident photon energy. By analyzing Ce 3$`d`$-edge RPES valence band spectra of CeFe2, CeNi2, and CeSi2 compounds systematically, we showed that this phenomenon originated from the Auger transition between different multiplet structures of $`\underset{¯}{3d}4f^2`$ electronic configuration in the intermediate state of the resonance process. Although we only showed the data on CeFe2, CeNi2, and CeSi2 in this paper, we found similar results for other Ce compounds as well, and expect this phenomenon to be universal for any Ce compounds. Hence, when the Ce $`f`$ spectral weights are obtained with the Ce $`3d`$-edge RPES technique, it is very important to avoid these overlapping Auger emissions contributing to the spectra. The most logical choice of the incident photon energy must be that corresponding to the lowest multiplet of the 3$`d`$-edge x-ray absorption spectra.

###### Acknowledgements.

One of us(EJCho) acknowledge the support by Korea Research Foundation Grant(KRF-2001-015-DP0178). This work is supported by the Korean Science and Engineering Foundation (KOSEF) through the Center for Strongly Correlated Materials Research (CSCMR) at Seoul National University. The research was performed under the support of Japan Synchrotron Radiation Research Institute(JASRI).

Figure Captions

Fig. 1 : x-ray absorption spectra of $`3d4f`$ transition region for CeFe2, CeNi2, and CeSi2. The symbols marked on the spectrum represent photon energies used for RPES shown in subsequent figures. The bar diagram at the bottom is the theoretical calculation result for the multiplet structure of $`\underset{¯}{3d}f^2`$ intermediate state taken from ref. .

Fig. 2 : Ce $`f`$ spectral weights near the Fermi level for CeFe2 obtained with various incident photon energies.

Fig. 3 : Difference spectra between the 4$`f`$ spectral weights at various photon energies and that at h$`\nu `$ = 881.76eV.

Fig. 4 : Ce 4$`f`$ spectral weights near the Fermi level for CeNi2 obtained at different photon energies. Shown in the bottom is the difference spectrum between the two.

Fig. 5 : Ce 4$`f`$ spectral weight near the fermi level for CeSi2 obtained at various photon energies. In the bottom the difference spectra are shown.


text:
# Bond-bending modes and stability of tetrahedral semiconductors under high pressure: a puzzle of AlN

## Abstract

Lattice vibrations of the würtzite-type AlN have been studied by Raman spectroscopy under high pressure up to the structural phase transition at $`20`$ GPa. We have shown that the widely debated bond-bending $`E_2^1`$ mode of w-AlN has an abnormal positive pressure shift up to the threshold of the phase transition, whereas in many tetrahedral semiconductors the bond-bending modes soften on compression. This finding disagrees with the results of ab initio calculations, which give a ”normal” negative pressure shift. Combination of high dynamical and low thermodynamical stability of AlN breaks the correlation between the mode Grüneisen parameters for the bond-bending modes and the transition pressure, which holds for CdS, InP, ZnO, ZnTe, ZnSe, ZnS, Ge, Si, GaP, GaN, SiC and BeO.

For the last four decades, a rather curious and at first sight surprising phenomenon of softening of tetrahedral semiconductors under compression has been reported in a number of experimental payne ; weinstein1 ; weinstein2 ; weinstein3 ; olego ; weinstein4 ; yakovenko ; perlin ; olijnik ; klotz  and theoretical dolling ; martin ; jex ; wendel ; yin ; chang ; nielsen ; chadi  works. Thus, Si, Ge, AIIIBV and AIIBVI semiconductors have pressure-sensitive phonon modes with negative pressure shifts weinstein2 ; yakovenko ; perlin ; klotz . These soft modes are shearing modes, involving bond bending in the first order of the strain yin ; chadi . Phonon frequency drop, more pronounced for high-Z materials, reaches about $`30`$% at the threshold of the pressure induced phase transitions, when covalent tetrahedral structures lose their stability and transform to more densely packed arrangements.

Earlier, the only known example of a bond-bending mode having a positive, although weak, pressure shift was being the ”pure bread” bond-bending $`E_2^1`$ mode of the low-Z würtzite-type BeO jephcoat . This result was being regarded as an exception until recent Raman measurements found a similar behavior for the $`E_2^1`$ mode of the würtzite-type AlN goni . Ab initio calculations goni ; gorczyca  did not catch this feature, giving a ”normal” negative pressure shift for the $`E_2^1`$ mode of AlN. This discrepancy deserves special attention, since the bond-bending elasticity is one of the most prominent manifestations of directional covalent bonding, and its pressure behavior should be well accounted for. The point to be made here is that AlN represents a complicated case of covalent versus ionic bonding karch : although its valence charge distribution is highly io nic gabe , AlN adopts the tetrahedrally coordinated würtzite structure under ambient conditions harrison . Aspiring to clarify the question, we have undertaken a complementary high-pressure Raman study of vibrational modes of the würtzite-type AlN up to its stability limit at about $`20`$ GPa. The pressure dependence of the low-frequency bond-bending $`E_2^1`$ mode is traced up to the threshold of the pressure induced phase transition for the first time.

The AlN samples were $`20`$ $`\mu `$m-thick crystals grown on the sapphire substrate by vapor phase epitaxy. Pressure was produced using the diamond-anvil pressure cell. Compressed helium and methanol-ethanol mixture were used as a pressure-transmitting medium in the first and in the second experimental run, respectively. Pressure was measured in situ by the ruby luminescence technique. The Raman spectra were measured using the THR-$`1000`$ triple spectrometer equipped with an OSMA detector (the first run), and the Dilor XY double spectrometer equipped with the CCD detector (the second run). An Ar+ laser ($`\lambda =514.5`$ nm) was used as a source of excitation. All spectra were recorded in the backscattering geometry at ambient temperature.

For the hexagonal würtzite structure with the space group $`P6_3`$mc ($`Z=2`$), a factor-group analysis predicts the following six sets of optical modes at the $`𝐤=0`$ fateley :

$$\mathrm{\Gamma }_{op}=A_1+2B_1+E_1+2E_2,$$

where $`A_1`$, $`E_1`$ and $`E_2`$ are Raman active modes, and $`B_1`$ modes are silent. $`A_1`$ and $`E_1`$ are also infrared active, and split into the longitudinal and transverse components (LO and TO). The lowest-frequency $`E_2^1`$ mode is a bond-bending mode.

The Raman spectrum of w-AlN has been measured previously under ambient conditions and analyzed in some detail, including the effects of polarization and anisotropy mcneil ; filippidis ; davydov . Our ambient pressure Raman frequencies $`\nu _{E_2^1}=249`$ cm-1, $`\nu _{A_1(TO)}=610`$ cm-1, $`\nu _{E_2^2}=657`$ cm-1, $`\nu _{E_1(TO)}=669`$ cm-1, $`\nu _{A_1(LO)}=890`$ cm-1, and $`\nu _{E_1(LO)}=910`$ cm-1 agree with very reliable data of Ref.mcneil ; filippidis ; davydov  within $`1\%`$. On increase in pressure, all Raman bands shift continuously to higher phonon energy with no intensity loss to about $`18`$ GPa. Above $`18`$ GPa, the bands weaken and disappear at about $`21`$ GPa in both experimental runs due to the phase transition to the rock salt structure ueno ; xia . Representative Raman spectra of w-AlN in the low-energy region as a function of pressure are shown in Fig. 1.

Fig. 2 compares measured and calculated pressure dependences of the Raman frequency for the $`E_2^1`$ mode. The mode pressure coefficients $`\nu _i^{}`$ calculated using the linear least-squares fit $`\nu _i=\nu _i^0+\nu _i^{}P`$, where $`\nu _i`$ is the frequency of the mode $`i`$ at pressure $`P`$, are listed in Table I comment1 . The pressure dependence of the $`E_2^1`$ frequency, measured in our experiment, is weak but apparently positive and linear at pressures to the threshold of the phase transition comment4 . This finding confirms the results obtained in goni , although our pressure slope of the $`E_2^1`$ frequency is somewhat lower (see Fig. 2). So, the positive pressure shift of the $`E_2^1`$ frequency of AlN should be regarded as a reliable experimental fact.

As we have mentioned above, ab initio calculations goni ; gorczyca  do not reproduce this feature of AlN, although the agreement between measured and calculated pressure shifts for the high-frequency bond-stretching $`E_2^2`$, TO and LO modes is very good (see Table I). Certainly, the value of the discussed effect for the $`E_2^1`$ mode is quite small, and the lack of ab initio calculations is not so crucial. Nevertheless, this lack is very instructive to realize a subtle pressure-sensitive balance between stabilizing and destabilizing contributions to the bond-bending force constant wendel ; yin ; martin  and may serve a further improvement of the theory.

Table II compiles the commonly used mode-Grüneisen parameters $`\gamma _i=d\mathrm{ln}\nu _i/d\mathrm{ln}V`$ obtained for the bond-bending modes $`i`$ in a series of tetrahedral compounds comment2 . Negative value for $`\gamma _i`$ is observed in each case, except for SiC, BeO and AlN. Although SiC has an essentially zero $`\gamma _i`$, its quadratic pressure coefficient is negative yakovenko . Thus, BeO and AlN appear to be the most stable materials with respect to the bond-bending mode on compression.

This bond-bending stability of AlN looks rather surprising, since its würtzite phase has a rather limited stability range up to $`21`$ GPa. Indeed, experimentally it has been found that the stability of tetrahedral structures with respect to the bond-bending modes correlates with their absolute stability under pressure in such a way that the frequency drop is faster for less stable compounds. Weinstein weinstein3 ; weinstein4  has discovered that for six diamond and zinc-blende structure semiconductors ZnTe, Ge, Si, ZnSe, ZnS and GaP there is a remarkable linearity between the mode Grüneisen parameter $`\gamma _{TA(X)}`$ for the bond-bending TA(X) mode and the transition pressure $`P_{tr}`$ for these materials. Ab initio calculations by Yin and Cohen yin  have given some insight into the nature of this correlation. They have shown that those individual contributions to the total crystal energy of Si and Ge, which stabilize the bond-bending TA(X) mode, are just the same ones, which preserve the diamond structure from rearrangement to more densely packed structures comment3 . The balance between stabilizing and destabilizing contributions drastically depends on the specific volume, destabilizing contributions becoming stronger under compression.

At present it is possible, following Weinstein and Zallen weinstein3 ; weinstein4 , to trace the $`\gamma _iP_{tr}`$ correlation to much higher pressures, using recent data for w-GaN perlin , SiC-$`6H`$ yakovenko ; yoshida  and w-BeO mori . The mode Grüneisen parameters determined for the bond-bending TA(L) phonons of InP, GaAs, ZnTe, ZnS, Ge, Si, GaP and for the bond-bending $`E_2`$ phonons of CdS, ZnO, AlN, GaN, SiC and BeO are plotted in Fig. 3 as a function of the transition pressure $`P_{tr}`$ ranging to $`140`$ GPa.

Fig. 3 shows that the $`\gamma _i`$ versus $`P_{tr}`$ points form a smooth curve deviating from linearity at $`\mathrm{P}_{tr}>20`$ GPa and approaching to zero at $`P_{tr}100`$ GPa. On the basis of this curve one may reasonably expect that materials with $`\gamma _i>0`$ should be at least as stable under pressure as SiC does. However, the point for AlN falls far afield from the common curve and tends to continue a linear dependence drawn by high-Z materials.

So, we see that a simple universal correlation between the pretransitional behavior of the bond-bending modes and the stability of tetrahedral semiconductors under high pressure does not exist. To our knowledge, the question why AlN combines a high dynamical and a low thermodynamical stability remains still open, and a detailed microscopic treatment, possibly along the direction indicated by Yin and Cohen yin , is required.

The authors wish to thank A. Dobrynin for growing the AlN crystals. E. V. Iakovenko is grateful to A. F. Goncharov for his assistance in Raman measurements.


text:
# Epitaxy of Fe3O4 on Si(001) by pulsed laser deposition using a TiN/MgO buffer layer

## I Introduction

Oxide thin films offer a large variety of physical properties useful for future electronic devices. Much attention has been paid to the ferroelectric, dielectric, and optical properties of perovskite materials with respect to non-volatile memory transistor applications or nonlinear optical devices Yu:00 . In the area of so-called spintronics, ferromagnetic oxides with high spin polarisation $`P`$ of the charge carriers at room temperature are desired for spin injection semiconductor devices as for example described by Datta and Das Datta:90 . While it has been shown that the spin lifetime in semiconductors is high as compared to metals Kikkawa:99 , one of the largest problems is to create a spin polarized population of the charge carriers in the semiconductor by injection of spin polarized currents. The use of diluted magnetic semiconductors - if existent \- is possible but limited due to the low Curie temperatures of the magnetic semiconductors so far Fiederling:99 ; Ohno:99 . Therefore, half-metallic materials with a full spin polarization of the charge carriers at the Fermi level and Curie temperatures well above room temperature are under consideration for spintronic devices. There are several candidates for materials with large $`P`$ close to 100% such as the Mn-based Heusler alloys, the oxide ferromagnets such as Fe3O4 Gupta:99  or CrO2, the doped manganites, and also the double perovskites Philipp:01 . We note that the use of these materials for ferromagnet-semiconductor hybrid devices requires the epitaxial growth of these materials on semiconductors. This is a great challenge with respect to materials technology. We further note that the resistivity mismatch between ferromagnetic metals and the semiconductors puts fundamental limits on the degree of spin polarization achievable with spin injection Schmidt:00 . The efficiency of spin injection, however, can be largely increased by introducing a tunnel barrier between the ferromagnetic metal and the semiconductor Rashba:00 ; Fert:01 .

Here, we report on the fabrication and structural characterization of magnetite (Fe3O4) epitaxial thin films on Si(001) substrates using a double buffer layer of TiN and MgO. MgO is also a candidate for an insulating barrier required for future tunneling injection devices. An important step is the in-situ substrate treatment by an infrared laser beam heating system. This allows for the controlled removal of the amorphous Si oxide surface layer prior to film deposition monitored by reflection high energy electron diffraction (RHEED). After substrate cleaning, the thin films are grown by ultra high vacuum (UHV) laser molecular beam epitaxy (L-MBE) Gross:2000a ; Klein:2001a . Also, the complete growth process is monitored by RHEED. We discuss in detail the growth of the buffer layer system consisting of TiN and MgO matching the magnetite lattice constants, and the epitaxial growth of Fe3O4 on top of the buffer layer system. The surface of the completed Fe3O4 thin films as well as the developing heterostructures is characterized after each process step by in situ atomic force microscopy (AFM). The structural properties of the completed hybrid heterostructures have been investigated by high-resolution x-ray diffraction (HRXRD).

## II Experimental Techniques

The heteroepitaxial thin film structures have been grown using laser molecular beam epitaxy (L-MBE) Gross:2000a ; Klein:2001a . The thin film UHV system consists of three separate sub-chambers for pulsed laser deposition, metallization, surface analysis and a load-lock. A sketch of the sub-chamber for the pulsed laser deposition (PLD) is shown in Fig. 1. A technical description of the systems has been given recently Gross:2000a ; Klein:2001a . Therefore, we only discuss those parts of the system relevant for the work reported here. For the present study, the two most important parts are (i) the (high-pressure) RHEED system and (ii) the infrared laser heating system used for the substrate heating.

### II.1 High pressure RHEED system

Many oxide materials have to be grown in a large oxygen partial pressure perturbing the RHEED analysis during growth due to the strong scattering of the high energy electrons from the oxygen molecules. Therefore, a special high pressure RHEED system has to be used, where the electrons travel only a short distance in the high pressure deposition atmosphere, and the rest of the beam path is pumped differentially Rijnders:97 . In our system a double differential pumping unit is used allowing for the operation of the RHEED system up to pressures of several 10 Pa (for details see Ref. Klein:2001a ).

### II.2 Laser heating of Si(001) substrates

An infrared laser heating system to heat substrates in a thin film deposition system to temperature above $`1200^{}`$C has first been introduced by Ohashi et al. Ohashi:99 . Such a system has several advantages with respect to thin film epitaxy. First, inside the vacuum chamber of the deposition system the substrate is heated directly by the infrared laser positioned outside the vacuum system. Therefore, the substrate heating is more efficient than in the case of a radiation heater positioned inside the vacuum chamber. In particular, the use of a radiation heater results in parts with a temperature well above the substrate temperature located inside the vacuum system. This in turn usually results in reduced vacuum conditions. In our case, using the infrared laser heating a substrate temperature $`T_S`$ well above 1200C can be reached at a vacuum in the 10-8 mbar range. Compared to the radiation heater setup in the same chamber, an improvement of the background pressure of about an order of magnitude is achieved. Second, the laser heating system allows for rapid changes of the substrate temperature, since only the substrate is heated, which has a much smaller heat capacity than the system radiation heater plus substrate. We realized changes of up to 500 K/min, only limited by the mechanical stability of the used substrate. Third, high temperatures up to the melting point of Si at 1410C can be reached without using an extensive amount of radiation shielding preventing the optical access to the sample surface.

In our setup the beam of a diode laser with a wavelength of 950 nm and a maximum power of 100 W is fed through an optical viewport and directly irradiates the backside of the Si(001) substrate. Using transparent substrates like SrTiO3, one has to take into account the temperature dependence of the absorption coefficient $`\eta `$. Here, a small amount of silver paint on the backside of the substrate is sufficient to absorb most of the laser power. This allows for the increase of the substrate temperature to values, where $`\eta `$ increases and thus enables substrate temperatures up to 1200C. In contrast, Si has high infrared absorption already at room temperature and therefore does not require any backside treatment. For temperature measurement a high resolution pyrometer working at the wavelength of $`1.5\mu `$m is used. The measured spot size is only 1 mm2, so even temperature inhomogeneities on the small substrates can be detected. Changes in the emission coefficient $`\epsilon `$ of the substrates at the measuring wavelength give an error in true temperature. Calibration measurements in an oven showed that this error with a big change in $`\epsilon `$ of 0.25 (worst case scenario) is only about $`\pm 3`$ K at $`400^{}`$C and $`\pm 20`$ K at $`1000^{}`$C. We note, however, that the for the deposition processes important reproducibility of the deposition temperature is with 5 K on the used substrates much better. Compared to a thermocouple the absolute exactness of this pyrometric measurement is good and the reproducibility strongly improved.

In the following we discuss the cleaning process for the surface of the Si(001) substrates. In our work we used commercial Si(001) single crystal wafers, which are one side polished and cut into 5 mm by 5 mm pieces. These pieces have been precleaned in aceton and isopropanol using an ultrasonic bath. However, it is unavoidable that a thin amorphous silicon oxide layer forms on the substrate surface. For the epitaxial growth of a thin film structure on the Si substrate, this amorphous oxide layer has to be removed. Here, an in situ cleaning process is of great advantage. For this purpose, the Si substrates have been laser heated in vacuum to $`T_S=900^{}`$C with a rate of 100 K/min and subsequently up to 1150C at 50 K/min. Then, the substrate was held for 10 min at this temperature. After that is was cooled down with a rate of 200 K/min to the deposition temperature of the first buffer layer (deposition of TiN at $`T_S=600^{}`$).

The whole high temperature treatment of the substrate surface was monitored by RHEED. The corresponding RHEED patterns are shown in Fig. 2. At low temperatures (Fig. 2a, $`T_S=600^{}`$), a large homogeneous background signal is visible due to diffuse scattering from the amorphous silicon oxide surface layer. Increasing the temperature, the signal from the (2$`\times `$1) surface reconstruction becomes visible within a lightly streaky pattern at about $`T_S=900^{}`$ (Fig. 2b). The surface reconstruction takes place only after oxygen and other impurities such as SiC have desorbed from the surface. It is interesting to note, that at temperatures above 1100C the (2$`\times `$1) surface reconstruction disappears. A possible explanation is that the dimer bonds of the (2$`\times `$1) surface displacement structure are dissolved at that temperature. Note that the diffuse background signal has disappeared, showing that the oxide has been completely removed from the substrate surface. After cooling down to $`1000^{}`$C, a stable and clean (2$`\times `$1) surface is obtained with strong diffraction spots on the Laue circle of 0th order, faint streaks, and Kikuchi lines. In summary, the process described above yields a clean Si(001) (2$`\times `$1) surface, as indicated by the presence of the additional diffraction spots, without any chemical surface treatment.

## III Experimental Results and Discussion

### III.1 Epitaxial growth of TiN, MgO, and Fe3O4

For the pulsed laser deposition of the TiN, MgO, and Fe3O4 thin films we used a 248 nm KrF excimer laser with an energy density at the target of 2-5 J/cm2 per shot, and a pulse repetition rate of 2-10 Hz, depending on the deposited material. Stoichiometric targets have been used for the growth of all films.

We first discuss the TiN films deposited in an Ar atmosphere of $`2.510^3`$ mbar and a pulse repetition rate of 5 Hz. As a function of substrate temperature we have found two different growth modes of the TiN(001) films on the Si(001) surface that have been cleaned prior to deposition as described above. As demonstrated by the optical micrographs shown in Fig. 3, above about $`T_S=650^{}`$C a three-dimensional Volmer-Weber (island) growth mode is obtained, whereas below this temperature a two-dimensional layer-by-layer growth mode is observed. For the Volmer-Weber growth mode the islands are oriented along the  direction of the Si(001) substrate. The size of the mostly quadratic islands was about $`3\mu \text{m}\times 3\mu \text{m}`$. It is not obvious why in the temperature range above $`650^{}`$C TiN grows with this preferential direction. We believe, that this is due to the 45 rotation of the dimer rows in the $`2\times 1`$ surface reconstruction. With the lattice constant of Si (diamond structure, $`a=0.543`$ nm) the distance of the Si atoms along the  direction is $`a/\sqrt{2}=0.384`$ nm resulting in a considerable lattice mismatch of 9.4% to TiN (cubic, $`a=0.424`$ nm). The transition into the island growth mode with increasing temperature itself is kinetic energy driven. For lower temperatures, a smooth two-dimensional growth mode is obtained. It is most likely that the growth occurs in the so-called 5-on-4-cube-on-cube bulk superstructure with a mismatch of about -2.4% Chowdury:94 ; Willmott:98 . In this superstructure 5 TiN unit cells overlap 4 Si(100) unit cells.

The RHEED pattern recorded during the growth of the TiN films is shown in Fig. 4. During the first 70 to 140 laser pulses corresponding to a film thickness of about 0.4 nm to 0.8 nm, TiN grows in a smooth superstructure resulting in a rich diffraction pattern (see Fig. 4a). With increasing film thickness the fine structure in the RHEED pattern disappears and clear diffraction spots are observed at positions expected for small island transmission (see Fig. 4b). This is consistent with an in situ AFM analysis, which yields an average surface roughness below 1 nm on a 1x1$`\mu `$m area at a film thickness of a few 10 nm.

After the growth of the TiN buffer layer the substrate temperature is lowered to $`T_S=330^{}`$C. Then, as a second buffer layer MgO is deposited at the same Ar pressure and a laser pulse repetition rate of 5 Hz. During the MgO growth the RHEED pattern becomes more streaky. This indicates a smoothening of the surface, which is consistent with in situ AFM measurements of the MgO surface. On a $`1\times 1\mu `$m2 area of the MgO surface the average roughness was about 0.6 nm. This shows that MgO is well suited as a buffer layer in the described system since it provides smooth interfaces, and therefore can also serve as insulator in a tunneling injection device.

On the completed TiN/MgO buffer layer system, an epitaxial magnetite film has been grown at $`T_S=330^{}`$C also in the same Ar atmosphere of $`2.510^3`$ mbar, at a laser pulse repetition rate of 2 Hz. These parameters have been found to result in the best RHEED intensity oscillations due to a layer-by-layer growth of Fe3O4 on MgO substrates Reisinger:03pre . Note that the RHEED pattern recorded during the growth of the magnetite film (see Fig. 5) is streaky but shows no indication of any transmission spots suggesting a smooth surface in contrast to the TiN layer. The small roughness of the surface of the magnetite film has been demonstrated directly by in situ AFM measurements performed after deposition. For an area of $`1\times 1\mu \text{m}^2`$ a root mean square (rms) roughness of only 0.3 nm has been found. The stripes in the RHEED pattern obtained for magnetite have only about half the spacing of the diffraction spots obtained for TiN and MgO. This is caused by the fact that the lattice constant of magnetite is about twice the one of TiN and MgO (see Figs. 4 and 5).

During the growth of magnetite on the TiN/MgO buffer layer system RHEED oscillations could be observed as shown in Fig. 6. From the total number of RHEED oscillations, the number of laser pulses per single oscillation and the film thickness determined from x-ray reflectometry (see next section), we find that four RHEED oscillations are obtained during the growth of a single unit cell thick layer. That is, the magnetite unit cell grows in four sub-unit-cell blocks as discussed in detail elsewhere Reisinger:03pre . After the growth of several unit cells the amplitude of the RHEED oscillations decreases and finally vanishes due to surface roughening. However, after an annealing step at the deposition temperature, the RHEED intensity recovers to about the same value observed before starting the deposition process. A key result is that a layer-by-layer growth mode of Fe3O4 can be realized on the TiN/MgO buffer layer system on Si(001) similar to the growth of Fe3O4 on MgO single crystalline substrates.

### III.2 Structural analysis by HRXRD

The structural properties of the thin films and the quality of epitaxy have been characterized using a four circle x-ray diffractometer. We have used an x-ray mirror and an asymmetric fourfold monochromator. This is required to have sufficient intensity and resolution in order to be able to distinguish between the diffraction peaks of the different layers which are expected to be close together.

The upper panel of Fig. 7 shows a $`\theta 2\theta `$ scan of the completed TiN/MgO/Fe3O4 multilayer system grown on Si(001). All layers have been grown in the ($`00\mathrm{}`$) orientation of the substrate yielding only ($`00\mathrm{}`$) peaks. Within the resolution of our x-ray diffraction system, no peaks of other chemical phases could be detected. The detailed positions of the ($`00\mathrm{}`$) peaks of the thin TiN and MgO buffer layers in the heterostructure are difficult to be determined (see lower panel of Fig. 7). Therefore, individual TiN films and TiN films with thicker MgO top layers have been grown for x-ray analysis. From these samples the following values for the $`c`$-axis parameter have been found: Fe3O4: $`c=0.847`$ nm, MgO: $`c=0.419`$ nm, and TiN: $`c=0.420`$ nm. The FWHM of the rocking curve of the magnetite (004) peak was about 0.7. This value is about ten times larger than the value measured for Fe3O4 films grown directly on MgO substrates under similar deposition conditions. This large difference can be attributed to the large lattice mismatch between Fe3O4 (as well as the TiN/MgO buffer layer) and Si compared to Fe3O4 and MgO (lattice mismatch of Fe3O4: -2.4% on Si(001) for the 5-on-4-cube-on-cube growth mode and only -0.31% on MgO). The large lattice mismatch results in strain relaxation effects and the formation of misfit dislocations increasing the mosaic spread of the films.

It is important to consider also the in-plane lattice constants in order to estimate whether or not the strain due to the lattice mismatch between the buffer layers and the substrate is (at least partially) relaxed. This can be done by several methods. A 2$`\theta \chi `$ mapping of the TiN (111) peak is shown in Fig. 8. Having determined the $`c`$-axis lattice parameter ($`c=0.420`$ nm) from the $`\theta 2\theta `$ scans, the in-plane lattice constant $`a`$ can be calculated for the tetragonal lattice of TiN as

$$a=\sqrt{\frac{h^2+k^2}{\left({\displaystyle \frac{2\mathrm{sin}\theta }{\lambda }}\right)^2\left({\displaystyle \frac{\mathrm{}}{c}}\right)^2}}.$$
(1)
Here, $`h`$, $`k`$ and $`\mathrm{}`$ are the Miller indices of the respective diffraction spot and $`\lambda `$ the wavelength of the Cu-Kα1 x-ray radiation used for the experiment. With $`c=0.420`$ nm, $`h=k=\mathrm{}=1`$, $`\theta =36.97^{}`$ and $`\lambda =0.15418`$ nm we obtain $`a0.422`$ nm. This value is close to the bulk value of TiN ($`a=0.424`$ nm). That is, most part of the TiN layer contributing to the intensity of the (111) peak is relaxed. From the mapping shown in Fig. 8 the in-plane lattice constant and the tetragonal distortion can be determined both from the $`\theta `$-value using eq.(1) and the $`\chi `$-value, which is given by the relation between the in-plane ($`a=b`$) and out-of-plane ($`c`$) lattice constants as

$$\frac{c\sqrt{2}}{a}=\mathrm{tan}\chi .$$
(2)
This relation directly follows from the symmetry of the tetragonal unit cell. For a nondistorted lattice (cubic symmetry), the $`\chi `$-value should be 54.74. As shown by Fig. 8, the measured value is 54.64 resulting in a $`c`$\- axis parameter of 0.42 nm, which is consistent with the value determined from the $`\theta 2\theta `$ scan. It is evident from Fig. 8 that the (111) peak is much broader in the $`\chi `$ than in $`2\theta `$ direction. This is caused by the measuring geometry using a line profile of the x-ray beam. Analyzing the magnetite (113) peak, we have determined the in-plane lattice constants of magnetite to $`a=b=0.832`$ nm by the same procedure. This value is close to the bulk value of $`a=0.835`$ nm. We note that the relaxation of the TiN buffer layer allows the growth of an almost strain free magnetite layer on Si(001).

We also performed $`\phi `$ scans to check the epitaxial relations between the different layers and the substrate. We found that the (111) TiN, (113) magnetite, and (111) Si peak in a $`\phi `$-scan (not shown) always appear at the same in-plane angle $`\phi `$ with the same fourfold symmetry. This shows that the in-plane lattice vectors of all layers are strictly aligned parallel to the Si substrate underlining the epitaxial nature of the obtained growth. We obtain the following epitaxial relations: $`_{Fe_3O_4}`$ $``$ MgO $``$ TiN $``$ Si.

The thickness of the individual layers within the TiN/MgO/Fe3O4 multilayer structure was determined by low angle x-ray reflectometry. The experimental results have been fitted using a simulation software based on dynamical scattering theory refsim . Fig. 9 shows a reflectometry measurement of a complete multilayer system consisting of a TiN/MgO/Fe3O4 multilayer on Si(001). Note the convincing agreement between the simulation and the measurement of the complicated heterostructure. From the refinement the layer thickness in this example is obtained to 43 nm for magnetite, 3.9 nm for MgO, and 12 nm for TiN. By fitting the data also the roughness of the surface as well as of the interfaces between the layers and the substrate could be estimated. The derived roughness values are in the range bewteen 0.3 and 0.7 nm. This is fully consistent with the roughness values obtained from the AFM anaylsis.

## IV Summary

Using pulsed laser deposition we have successfully grown magnetite epitaxial thin films on Si(001) substrates using a double buffer layer system consisting of TiN and MgO. For the reproducible cleaning of the Si substrate we have used an in situ laser heating process with a maximum temperature of 1150C. The substrate treatment and the whole growth process was monitored directly by RHEED allowing the determination of the surface structure at each process step. By x-ray diffraction and AFM, the good crystalline quality of the multilayers and the small roughness of the interfaces were verified. Due to their high Curie temperature and predicted large spin polarization magnetite thin films on Si substrates are highly interesting for potential applications in spintronics, in particular for spin injection into semiconductors. The compatibility of the presented multilayer structure with the insulating material MgO may provide a way to fabricate a suitable tunnel barrier to overcome the resistivity mismatch between ferromagnetic metals and semiconductors.

This work was supported in part by the Deutsche Forschungsgemeinschaft (project No.: Al/560) and the BMBF (project No.: 13N8279).