text:
# Starbursts from strong compression of galactic molecular clouds due to the high pressure of the intracluster medium

## 1 Introduction

Recent spectrophotometric and morphological studies of cluster galaxies by large ground–based telescopes and the $`Hubble`$ $`Space`$ $`Telescope`$ ($`HST`$) have revealed that the star formation histories are very different between different morphological and luminosity classes and can be closely associated with galactic morphological evolution driven by cluster-related physical processes (e.g., Lavery & Henry 1994; Couch et al. 1994, 1998; Abraham et al. 1996; Balogh et al 1999; Dressler et al. 1999). For example, Couch et al. (2001) revealed that star formation of disk galaxies in distant clusters are globally and uniformly suppressed. Poggianti et al. (1999) found that galaxies with “e(b)” spectral type, which are the most likely to be starbursting galaxies, generally belong to the low-luminosity dwarf populations in clusters. It has been a longstanding and remarkable problem what mechanisms are responsible for the observed enhancement/decline of galactic star formation rates in clusters.

The dynamical and hydrodynamical effects of the hot, high-temperature ICM have long been suggested to play decisive roles, not only in transforming galactic morphological properties but also in controlling star formation rates in galaxies within clusters (e.g., Gunn & Gott 1972; Gavazzi & Jafffe 1985; Bothun& Dressler 1986; Evrard 1991). Fujita & Nagashima (1999) demonstrated that star formation rates and photometric properties of galaxies in clusters can be changed as a result of the ram pressure effects on galactic molecular clouds. Using three dimensional SPH/N-body simulations, Abadi et al. (1999) found that although the ram pressure of the ICM can efficiently strip interstellar HI gas from spirals, such ram pressure stripping alone can not abruptly truncate star formation in spirals. However, these previous phenomenological/numerical models have not investigated the effects of the ICM on individual molecular clouds with internal structures, kinematics, and chemical abundance, even though star formation is ongoing within individual molecular clouds. Therefore, it is still unclear how the star formation rate and efficiency in individual galactic molecular clouds is influenced by the hot ICM (and thus how the galactic star formation can be controlled by the ICM).

The purpose of this Letter is to first demonstrate that the high gaseous pressure of the ICM significantly changes the internal structure of a self-gravitating molecular gas cloud in a cluster disk galaxy and consequently trigger a burst of star formation within the cloud. In particular, we investigate (1) the star formation caused within gas clouds in a disk galaxy due to the strong compression of the gas by the ICM, and (2) its dependence on the density and the temperature of the ICM and the velocity of the galaxy relative to the ICM. We show that both ram pressure and static pressure from the ICM on the galactic gas clouds are important in compressing the clouds and thus in triggering starbursts within the clouds. We emphasize that although the ICM can be responsible for the stripping of HI gas within disks and diffuse halo gas in cluster galaxies and thus for the truncation of star formation (e.g., Abadi et al. 1999; Bekki et al. 2002), it could also be closely associated with starburst activity in cluster/group environments. Finally, this proposed new mechanism for triggering starbursts is discussed in a variety of different contexts of cluster galaxy evolution.

## 2 Our Model

By using TREESPH code with star formation methods (Bekki 1997), we numerically investigate the hydrodynamical effects of the ICM on a self-gravitating molecular gas cloud orbiting within a spiral galaxy under the gravitational influence of the galaxy and its surrounding dark matter halo. The cloud is represented by 20,000 SPH particles and the initial cloud mass ($`M_{\mathrm{cl}}`$) and size ($`r_{\mathrm{cl}}`$) are set to be $`10^6M_{}`$ and 97 pc, respectively, which are consistent with the mass-size relation observed by Larson (1981). The cloud is assumed to have an isothermal radial density profile with $`\rho (r)1/(r+a)^2`$, where $`a`$ is the core radius of the cloud and set to be $`0.2r_{\mathrm{cl}}`$. An isothermal equation of state with a sound speed of $`c_\mathrm{s}`$ is used for the gas, and $`c_\mathrm{s}`$ is set to be 4 km s-1 (consistent with the prediction from the virial theorem) for models with $`M_{\mathrm{cl}}=10^6M_{}`$.

Each SPH particle is subject to the gravitational forces from the fixed spiral potential that is assumed to have three components: a dark matter halo, a disk, and a bulge. We assume a logarithmic dark matter halo potential ($`\mathrm{\Phi }_{\mathrm{halo}}=v_{\mathrm{halo}}^2\mathrm{ln}(r^2+d^2)`$), a Miyamoto-Nagai (1975) disk ($`\mathrm{\Phi }_{\mathrm{disk}}=GM_{\mathrm{disk}}/\sqrt{R^2+(a+\sqrt{z^2+b^2})^2}`$), and a spherical Hernquist (1990) bulge ($`\mathrm{\Phi }_{\mathrm{bulge}}=GM_{\mathrm{bulge}}/(r+c)`$), where $`r`$ is the distance from the center of the spiral, $`d`$ = 12 kpc, $`v_{\mathrm{halo}}`$ = 131.5 km $`\mathrm{s}^1`$, $`M_{\mathrm{disk}}`$ = $`10^{11}`$ $`M_{}`$, $`a`$ = 6.5 kpc, $`b`$ = 0.26 kpc, $`M_{\mathrm{bulge}}`$ = 3.4 $`\times `$ $`10^{10}`$ $`M_{}`$, and $`c`$ = 0.7 kpc. This set of parameters is reasonable and realistic for the Galaxy. The cloud is initially within the disk plane (corresponding to the $`x`$-$`y`$ plane) and located at ($`x`$, $`y`$, $`z`$) = ($`R_g`$, 0, 0). The initial velocity of the cloud is set to be ($`v_\mathrm{x}`$, $`v_\mathrm{y}`$, $`v_\mathrm{z}`$) = (0, $`v_\mathrm{c}`$, 0), where $`v_\mathrm{c}`$ is the circular velocity at the position. We show the results for $`R_g=8.5`$ kpc (corresponding to the solar neighborhood) in this paper, because we find that the star formation within clouds does not depend strongly on $`R_g`$.

Each SPH particle is also subject to the hydrodynamical force of the ICM whose strength depends on the parameters of the ICM. The ICM is represented by 8,824 SPH particles and and an isothermal equation of state with temperature, $`T`$ (or sound velocity of $`c_\mathrm{h}`$), and density, $`\rho `$. The cloud – which is in a spiral that orbits the cluster core with a velocity, $`V_{\mathrm{rel}}`$, relative to the ICM – is subject to both ram pressure ($`P_{\mathrm{ram}}`$ $``$ $`\rho V_{\mathrm{rel}}^{}{}_{}{}^{2}`$) and thermal static pressure ($`P_\mathrm{s}`$ $``$ $`\rho c_{\mathrm{h}}^{}{}_{}{}^{2}`$). The ICM particles are uniformly distributed in an annular ring with size $`R_\mathrm{g}`$, width $`6R_{\mathrm{cl}}`$, and thickness $`6R_{\mathrm{cl}}`$ (in the $`z`$ direction) on the disk. We give all ICM particles an initial velocity of (0, $`V_{\mathrm{rel}}`$, 0), so that we can represent the hydrodynamical effects that the ICM has on the spiral when it passes through the cluster core. Periodic boundary conditions are adopted only for these ICM particles. We here stress that for the adopted total number of SPH particles ($``$30000), the simulations can not fully resolve the contact surface between the cold GMC and the hot ICM (accordingly thermal conduction between these gas and the cloud evaporation can not be fully investigated). The star formation rate in our simulations can be therefore overestimated (as discussed later in this paper).

A gas particle in a cloud is converted into a collisionless stellar particle if (1) the local dynamical time scale \[corresponding to $`(4\pi G\rho _i)^{0.5}`$, where $`G`$ and $`\rho _i`$ are the gravitational constant and the density of the gas particle, respectively\] is shorter than the sound crossing time (corresponding to $`h_i/c_\mathrm{s}`$, where $`h_i`$ is the smoothing length of the gas), and (2) the gaseous flow is converging. This method thus mimics star formation due to Jeans instability in gas clouds. In the model without ICM effects, star formation does not occur at all in this star formation method. By changing the values of the three parameters (i.e., $`T`$, $`\rho `$, and $`V_{\mathrm{rel}}`$), we investigate the effects of the ICM on the cloud. We show the results for the models with $`0.01T_0`$ $``$ $`T`$ $``$ $`T_0`$, with $`0.1\rho _0`$ $``$ $`\rho `$ $``$ $`\rho _0`$, and $`0.1V_0`$ $``$ $`V_{\mathrm{rel}}`$ $``$ $`V_0`$, where $`T_0`$, $`\rho _0`$, and $`V_0`$ are 8 keV, 5.6 $`\times `$ $`10^{27}`$ g cm-3, and 1000 km s-1. These values of $`T_0`$ and $`\rho _0`$ correspond to the observed values of the ICM in the central region of the Coma cluster (Briel et al. 1992). We mainly show the result of the “fiducial model” with $`T`$ = $`T_0`$, $`\rho `$ = $`0.1\rho _0`$, and $`V_{\mathrm{rel}}`$ = 1000 km s-1, because this model shows the typical behavior of ICM-induced star formation in a gas cloud. The gas cloud can not be suddenly stripped because of the higher initial column density during simulations (See also Fujita & Nagashima 1999).

## 3 Results

Figures 1 and 2 describe how a burst of star formation within a gas cloud can be triggered by the external gas pressure of the ICM in the fiducial model. In this model with moderately strong ram pressure, the high pressure of the ICM can continue to strongly compress the cloud without losing a significant amount of gas from the cloud (i.e., only 3% of the initial mass can be stripped within 14 Myr). As the strong compression proceeds, the internal density/pressure of the cloud can rise significantly. However, the self-gravitational force of the cloud also becomes stronger because the cloud becomes progressively more compact during the ICM’s compression. Therefore, the internal gaseous pressure of the cloud alone becomes unable to support itself against the combined effect of the external pressure from the ICM and the stronger self-gravitational force. As a result of this, the cloud’s collapse initially induced by the high external pressure from the ICM can continue in a runaway manner.

Due to the rapid, dissipative collapse, the gaseous density of the cloud dramatically rises and consequently star formation begins in the central regions of the cloud. The star formation rate increase significantly from $`0.1`$ $`M_{}`$ yr-1 (for the first 5 Myr) to 0.6 $`M_{}`$ yr-1 (8 Myr after the start of the cloud’s collapse). This corresponds to global galactic star formation rate of $`3060`$ $`M_{}`$ yr-1 for galaxies with total H2 mass similar to that of the Galaxy. About 81% of the gas is converted into stars within 14 Myr to form a stellar system. Because of the “implosive” formation of stars from strongly compressed gas, the developed stellar system is strongly self-gravitating and compact. This result implies that high external pressure from the ICM is likely to trigger the formation of bound, compact star clusters rather than unbound, diffuse field stars.

Parameter dependences are found to be complicated and can be summarized as follows (see Figure 3): Firstly, for the high ICM temperature ($`T`$ = $`T_0`$), large relative velocity ($`V_{\mathrm{rel}}`$ = $`V_0`$) models, the star formation efficiency (hereafter SFE) within a cloud is higher for the model with the lower ICM density ($`\rho `$). This is because in the higher $`\rho `$ model, the ram pressure of ICM strips a larger amount of gas from the cloud and thus decreases the total amount of gas that can be converted into stars. Secondly, for the models with lower $`\rho `$ (= $`0.1\rho _0`$) and lower $`T`$ (= $`0.1T_0`$), the SFE is higher for the model with the larger $`V_{\mathrm{rel}}`$. This suggests that the moderately strong ram pressure of the ICM can play a role in enhancing the star formation rate in a gas cloud. Thirdly, for models in which ram pressure has a negligible effect on the gas cloud (i.e, high density, $`\rho `$ =$`\rho _0`$, and very low ICM relative velocity, $`V_{\mathrm{rel}}`$ = 0.01$`V_0`$), the SFE is likely to be higher for the model with higher $`T`$. This is because the higher static (thermal, external) pressure of the ICM in the model with higher $`T`$ can more strongly compress a cloud and thus produce regions of high gaseous density within. Fourthly, for models in which ram pressure is very strong ($`\rho `$ = $`\rho _0`$ and $`V_{\mathrm{rel}}`$ = $`V_0`$), there is no clear dependence on $`T`$.

## 4 Discussion and conclusions

By assuming that galactic molecular clouds can be exposed to the hot ICM, we have demonstrated that the high pressure from the ICM can trigger efficient star formation in gas clouds in cluster spirals. However, the diffuse outer halo gas and the disk HI gas of a cluster spiral could prevent the molecular gas clouds from being directly exposed to the hot ICM. Using semianalytic models of galaxy formation, Okamoto & Nagashima (2003) demonstrated that molecular gas clouds can be consumed up by star formation in cluster spirals before the ram pressure effects become significant for their evolution. Furthermore, evaporation of gas clouds by heat conduction (Cowie & McKee 1977) could reduce significantly the gas mass that can be converted into stars, because the time scale of cloud evaporation is an order of $`10^610^8`$ years for the adopted parameters of gas clouds in the present study. For the adopted cloud parameters, the fraction of gas that can be evaporated before star formation within the gas can be estimated as 12 %. Therefore, the proposed ICM-induced starbursts in molecular gas clouds may be more likely for galaxies which only have small amounts of HI and halo gas yet have enough molecular gas. In spite of this strong limitation in the applicability of the proposed mechanism to the starburst phenomenon within clusters, we suggest that the present results may still have the following three important implications for cluster galaxy evolution.

The first is that if the ICM triggers a starburst in a disk galaxy rich in molecular gas, actively star-forming regions can be widely spread throughout the disk. The expected distribution of starburst regions within a disk is in striking contrast to the nuclear starbursts suggested to be triggered by tidal interactions/mergers between galaxies and global cluster tidal fields (e.g., Byrd & Valtonen 1990). Furthermore, the starburst mechanism proposed here for cluster disk galaxies differs from that for mergers in that the thin disk component remains intact after the starburst in the former whereas it is destroyed to form a thick disk in the latter (e.g., Bekki 1998). Therefore, widely spread starburst (or poststarburst) activity within the thin disk component of a cluster disk galaxy is the signature to look for in identifying this ICM-gas cloud interaction. Spectroscopic observations that can reveal the detailed distributions of starburst/poststarburst regions within a disk (e.g., the GMOS Integral Field Unit on the 8m Gemini telescope) will reveal such ICM-induced starburst/poststarburst populations in distant clusters.

The second is that a larger fraction of starburst/poststarburst populations are expected in merging clusters, where a higher ICM temperature (and pressure) is predicted (e.g., Roettiger et al. 1997). However, it is unclear (both observationally and theoretically) whether cluster merging triggers starbursts in cluster disk galaxies and what mechanisms are responsible for the formation of starburst galaxies in merging clusters. Although some observational evidence has been found for merging clusters having a larger number of starbursts/poststarburst candidates both at low (Caldwell et al. 1993; Caldwell & Rose 1997) and intermediate redshifts (e.g., Owen et al. 1999), Venturi et al. (2000) found no such clear evidence in the merging clusters Abell 2125 and 2645. Bekki (1999) demonstrated that the time-dependent tidal gravitational field of merging clusters can trigger starbursts in their galaxies whereas Fujita et al. (1999) concluded that ram pressure stripping can strongly suppress the star formation of galaxies in merging clusters. More extensive observational investigations of (1) a correlation between the fractional content of starburst/poststarburst galaxies and the incidence/non-incidence of cluster merging, and (2) the spatial distribution of starburst/poststarburst galaxies within merging clusters, will address the above two unresolved problems in a more quantitative way.

The third is that the origin of isolated intracluster compact HII regions recently discovered in the Virgo cluster (e.g., Yoshida et al. 2002; Arnaboldi et al. 2003; Gerhard et al. 2003) can be closely associated with ICM-induced efficient star formation in gas clouds drifting in the cluster. The ICM-induced efficient star formation might be highly likely to be seen in the intracluster molecular gas clouds stripped from cluster spirals through interactions between them (and between them and the cluster global tidal fields), because these intracluster gas clouds are directly exposed to ICM. We suggest that the total number and spatial distribution of isolated intracluster HII regions in the Virgo cluster provide valuable information on the cluster’s dynamical history of galaxy interaction/merging that have been responsible for stripping of molecular clouds within disks. Our results also imply that the formation of intracluster globular clusters (West et al. 1995) and planetary nebulae (Arnaboldi et al 1996) could result from the past interaction between stripped gas clouds and ICM.

The present results imply that some fraction of molecular clouds in cluster spirals could disappear from their disks, not by tidal or ram pressure stripping but by rapid consumption due to ICM-induced star formation. They also imply that the initial mass function (IMF) of a gas cloud with ICM-induced star formation could be different from the normal Salpeter-like IMF, because the ICM-induced star formation can occur preferentially in the central regions of gas clouds, where more massive stars are likely to be formed (e.g., Murray et al. 1996). Our future, more sophisticated simulations which will include chemical evolution, magnetic fields, dynamical evolution of hierarchical/fractal structures within a cloud, and feedback effects from massive stars and supernovae will address the total amount of molecular gas consumed by ICM-induced star formation and its IMF in a more quantitative way.

We are grateful to the anonymous referee for valuable comments, which contribute to improve the present paper. K.B. and W.J.C. acknowledge the financial support of an Australian Research Council Discovery Program grant throughout the course of this work.


format: text

text:
# Sierpinski signal generates 1/𝑓^𝛼 spectra

## Definition of the model:

The dynamics of the so-called Sierpinski automaton is related to the generation law of the Pascal triangle. This pattern can be generated by the following simple one-dimensional cellular automaton: We consider a linear array of sites (or spins) $`x_i(t)`$ which can take the values 0 or 1 at discrete time steps $`t`$. We restrict ourselves to the special initial condition, that for $`t=0`$ only one spin is different from all others:

$`x_0(0)=1\mathrm{and}_{i0}x_i(0)=0.`$ (1)
The dynamics is defined by the following next-neighbor interaction:

$`x_i(t+1)=[x_{i1}(t)+x_{i+1}(t)]\mathrm{mod}2,`$ (2)
$`i[\mathrm{},\mathrm{}]`$, $`x_i(t)\{0,1\}`$ at discrete time $`t`$. In the context of catalytic processes otterstedt98 ; dress84plath ; liauw96  a simplified chemical interpretation of this rule reads: A catalytic process is stopped when too less (i.e. no) or too many (i.e. 2) neighbor sites are active. A catalytic process is initiated (or continued) when only one neighbor site is active. This can origin from a minimal catalysis temperature combined with a local self-limiting reaction rate. The spatiotemporal evolution is obtained from the rows of the Pascal triangle by applying modulo 2 elementwise,

$`\begin{array}{ccccccccccccccccccc}& & & & 1& & & & & & & & & & 1& & & & \\ & & & 1& & 1& & & & & & & & 1& & 1& & & \\ & & 1& & 2& & 1& & & \stackrel{\text{modulo 2}}{}& & & 1& & 0& & 1& & \\ & 1& & 3& & 3& & 1& & & & 1& & 1& & 1& & 1& \\ 1& & 4& & 6& & 4& & 1& & 1& & 0& & 0& & 0& & 1\end{array}`$ (8)
which is also known as the Sierpinski gasket. The Sierpinski gasket is a well-known self-similar structure (with point dimension $`\mathrm{ln}3/\mathrm{ln}2`$ in $`x`$ for $`t\mathrm{}`$) that is obtained by twofold replication of the first four rows to the subsequent four rows, and iteration of this process with the whole triangle.

Instead of considering the fractal pattern itself, we look at a scalar observable that can be compared to experimental time series. Before considering the spectrum, we briefly sketch a direct solution kummer1852  and illustrate the analogy to a formal language approach. The row sum (or total (in-)activity) over space at time $`t`$, defined by

$`X(t):={\displaystyle \underset{i}{}}x_i(t),`$ (9)
is referred to as Sierpinski signal. The Sierpinski automaton rule then generates a time series $`X(t)`$ (fig. 1) starting from $`t=0`$ with

$`1,2,2,4,2,4,4,8,2,4,4,8,4,8,8,16,\mathrm{}`$ (10)
Being interested in an analytic expression for (10), we first note that $`X(t)`$ can be generated up to $`t=2^N1`$ from the start sequence $`u_0=(1)`$ by $`N`$ iterations of the sequence replication rule

$`u_nu_{n+1}=(u_n,2u_n).`$ (11)
Obviously $`X(t)`$ takes only powers of 2, so we consider $`(\mathrm{ln}X(t))/\mathrm{ln}2`$, which starts as

$`0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,\mathrm{}`$ (12)
Strikingly, this appears to be the number of ones in

$`0,01,10,11,100,101,110,111,1000,1001,`$
$`1010,1011,1100,1101,1111,\mathrm{},`$ (13)
i.e. the binary decomposition of the time variable $`t`$ starting with $`t=0`$. Therefore the observable is given by

$`X(t)=2^{\mathrm{CrossSum}\left(\mathrm{Binary}\left(t\right)\right)},`$ (14)
which is no longer recursively defined. This analytic solution has already been known in 1852 by Kummer kummer1852  in a number-theoretic context.

A convenient closed expression can be obtained by expressing time by a number of spins

$`t={\displaystyle \underset{j=0}{\overset{N1}{}}}\sigma _j2^j\sigma \{0,1\}.`$ (15)
and expressing $`X(\{\sigma _i\})`$ from the same configuration as

$`X(t)=X\left({\displaystyle \underset{j=0}{\overset{N1}{}}}\sigma _j2^j\right)=2^{_{j=0}^{N1}\sigma _j}.`$ (16)
By these expressions for $`t(\{\sigma _i\})`$ and $`X(\{\sigma _i\})`$, we have a parametric expression parameterized by a set of $`N`$ spins for all $`X(t)`$ with times up to $`t=2^N1`$.

## Spectrum of the Sierpinski signal:

The periodogram $`X(\omega )`$ of the time signal now is calculated analytically. The binary time decomposition allows a Fourier transformation of $`X(t)`$ fairly direct from the definition,

$`X(\omega )`$ $`=`$ $`{\displaystyle \underset{t=0}{\overset{2^N1}{}}}\mathrm{e}^{\mathrm{i}\omega t}X(t)={\displaystyle \underset{\sigma _0}{}}\mathrm{}{\displaystyle \underset{\sigma _{N1}}{}}\mathrm{e}^{\mathrm{i}\omega t(\{\sigma _i\})}X(t(\{\sigma _i\}))`$ (17)
$`=`$ $`{\displaystyle \underset{\sigma _0}{}}\mathrm{}{\displaystyle \underset{\sigma _{N1}}{}}{\displaystyle \underset{j=0}{\overset{N1}{}}}\mathrm{exp}\left(\sigma _j(\mathrm{i}\omega 2^j+\mathrm{ln}2)\right)`$
$`=`$ $`{\displaystyle \underset{j=0}{\overset{N1}{}}}{\displaystyle \underset{\sigma _j}{}}\mathrm{exp}\left(\sigma _j(\mathrm{i}\omega 2^j+\mathrm{ln}2)\right)`$
$`=`$ $`{\displaystyle \underset{j=0}{\overset{N1}{}}}\left(1+\mathrm{exp}(\mathrm{i}\omega 2^j+\mathrm{ln}2)\right),`$
where all sums over $`\sigma _j`$ are taken over the two possible values $`\sigma _j=0`$ and $`\sigma _j=1`$. We now calculate the periodogram’s power spectrum, i. e. $`|X(\omega )|^2`$. The absolute value of $`X(\omega )`$ simplifies to a trigonometric product which the logarithm converts into a sum,

$`\mathrm{ln}|X(\omega )|^2={\displaystyle \underset{j=0}{\overset{N1}{}}}\mathrm{ln}[5+4\mathrm{cos}(\omega 2^j)],`$ (18)
showing a rugged fine structure as shown in fig. 4. A rough estimate of the sum in eq. (18) is obtained approximating the sum by an integral ($`y:=\omega 2^j`$),

$`\mathrm{ln}|X(\omega )|^2`$ $``$ $`{\displaystyle _0^{N1}}\mathrm{ln}[5+4\mathrm{cos}(\omega 2^j)]dj`$ (19)
$`=`$ $`{\displaystyle \frac{1}{\mathrm{ln}2}}{\displaystyle _\omega ^{\omega 2^{N1}}}{\displaystyle \frac{\mathrm{ln}[5+4\mathrm{cos}y]}{y}}dy.`$ (20)
As $`\mathrm{ln}(5+4x)\mathrm{ln}(5)+\frac{4}{5}x`$ for $`|x|1`$, we obtain

$`\mathrm{ln}|X(\omega )|^2{\displaystyle \frac{\mathrm{ln}5}{\mathrm{ln}2}}{\displaystyle _\omega ^{\omega 2^{N1}}}{\displaystyle \frac{\mathrm{d}y}{y}}+{\displaystyle \frac{4}{5\mathrm{ln}2}}{\displaystyle _\omega ^{\omega 2^{N1}}}{\displaystyle \frac{\mathrm{cos}(y)}{y}}dy.`$
The integral over the integral cosine is nearly independent of the upper boundary for high values of the boundary. Thus, we can substitute the upper boundary $`\omega 2^{N1}`$ by some $`N`$-dependent constant, say $`c_N1`$. Finally, substituting the cosine by one yields immediately a rough approximation of eq. (18):

$`|X(\omega )|^2c_N^{}\omega ^{4/(5\mathrm{ln}2)}\omega ^{1.15}.`$ (21)
Due to the increase of the mean of $`X(t)`$, spectral estimation from the periodogram $`X(\omega )`$ (implying a periodical extension in time domain) has to be discussed carefully. As $`_{t=0}^{2^N1}X(t)=3^N`$, the average increase is $`X(t)_\mathrm{L}=X(t)_{\{0\mathrm{}2^N1\}}=t^\beta `$ with $`\beta =\mathrm{log}(3/2)/\mathrm{log}(2)0.585`$. Further, the signal exhibits an increasing variance $`X(t)^2_{\{0\mathrm{}2^N1\}}=t^\gamma `$ with $`\gamma =\mathrm{log}(5/2)/\mathrm{log}(2)1.32`$. Consequently, we investigate two variants of $`X(t)`$: A suitable per definitionem mean-free sum signal defined by $`Y(t)=X(t)(1+\beta )t^\beta `$, and a mean-free signal with non-increasing variance

$`Z(t)=Y(t)/Y(t)^2_{\{0\mathrm{}t\}}^{1/2}.`$ (22)
The spectrum of $`Y(\omega )`$ can be directly obtained from eq. (17) and the evaluation of the Fourier transform of $`t^\beta `$:

$`Y(\omega )=X(\omega )(1+\beta )(t^\beta )`$ (23)
The power spectrum of the periodically extended function $`t^\beta `$ decays (for small values of the frequency) like a power law with an exponent of approximately -2. Thus, the decay is much stronger than $`X(\omega )`$ and the power spectra of $`Y(t)`$ and $`X(t)`$ deviate only slightly. Hence, it follows that $`|Y(\omega )|^2|X(\omega )|^2\omega ^\alpha `$. Similarly, $`|X(\omega )|^2`$ also estimates $`|Z(\omega )|^2`$. We now compare these results with numerically applied discrete Fourier transformation. The power spectrum is fitted (least-squares) in fig. 4 by a power law with exponent $`\alpha `$ about 1.11, being in good agreement with the analytical result $`\alpha 1.15`$ of eq. (21).

If one measures a power spectrum experimentally, this may generically be done by observation of resonances, where the system is coupled with a tuneable oscillator of given frequency and finite bandwidth. Therefore it is quite natural to consider an averaged spectrum. The (incommensurable) averaging procedure applied in fig. 4 smoothes the peaks at $`\omega =2^k`$. The peak amplitudes decay as the average spectrum itself. For commensurable averaging (inset of fig. 4) the peaks at $`\omega =2^k`$ disappear completely.

Note that the spectra $`|X(\omega )|^2`$, and $`|Y(\omega )|^2`$, from FFT for $`T=2^{22}`$ (not shown) also display $`1/f^\alpha `$ behaviour, with exponents $`\alpha _X=1.12`$, and $`\alpha _Y=1.11`$ , respectively. Moreover, we have used the method of a sliding window that normalize the fluctuations Hu2001  of the detrended signal $`Y(t)`$, i.e. $`\stackrel{~}{Z}(t)=Y(t+l)/Y^2_{\{tl+1,t+l\}}^{1/2}`$. For different values of the window width $`2l`$, the power spectra exhibit power law behaviour with exponents of about $`\alpha 1.1`$. Thus $`1/f^\alpha `$ spectra appear to be a robust property of Sierpinski signals.

## Amplitude distribution:

Many systems exhibiting $`1/f`$ noise possess a Gaussian amplitude distribution voss1978 . In this paragraph we calculate the amplitude distribution $`H_N(2^k)`$ of $`X(t)`$ analytically where $`H_N(2^k)`$ denotes the frequency occurrence of $`X=2^k,k=0,1,\mathrm{}`$ for a signal length of $`T=2^N`$. The number of $`2^k`$s in the signal up to $`t=2^N1`$ is the number of $`2^k`$s plus the number of $`2^{k1}`$s in the signal sequence up to $`t=2^{N1}1`$, i.e. $`H_N(2^k)=H_{N1}(2^k)+H_{N1}(2^{k1})`$. For the boundary condition $`H_N(1)=H_1(2)=1`$ we obtain a sum of binomial coefficients $`H_N(2^k)=_{j=1}^N\left(\genfrac{}{}{0pt}{}{j}{k}\right)=\left(\genfrac{}{}{0pt}{}{N+1}{k+1}\right),k1`$, which simplifies for $`N1`$ to $`2^N`$ times the Gaussian distribution:

$`H_N(2^{k1}){\displaystyle \frac{2^N}{(\pi N/2)^{1/2}}}e^{\frac{(kN/2)^2}{N/2}}.`$ (24)
Hence, the amplitude distribution of the occurrence of powers of 2 is Gaussian for fixed $`N`$. The (averaged) amplitude distributions for $`Y(t)`$ and $`Z(t)`$ differ from $`H_N(2^k)`$ but possesses a similar shape as $`H_N(2^k)`$. Note that the variance distributions for $`X(t)`$, $`Y(t)`$, and $`Z(t)`$, are not Gaussian but well defined by eq. (24).

As a final point, numerical simulations show that the averaged signals are robust against noise, i.e. initial conditions with more than a single 1.

In analogous situations in less simply defined systems, power laws have also been observed in spatial spectra of the scum on fluid surfaces and in the random baker map AntonsenPRL95 , and in the temporal spectra in dissipative dynamics governed by the Lorenz equations pikovskyPRE95  being related to the Thue-Morse sequence pikovskyPRE95 ; zaksPRL96 ; zaksJSP97 ; zaksPRE01 . In fact the Thue-Morse dynamics $`1(1,1),1(1,1)`$ itself maps on the string replication rule $`u_n(u_n,(1)u_n)`$ for generation of the spatial sequence. Being not equivalent, but of striking similarity to eq. (11) for generation of the temporal Sierpinski signal series, the Sierpinski dynamics itself does not follow a replication rule. While the analytic solution of the Thue-Morse sequence is $`(1)^{\mathrm{CrossSum}(\mathrm{Binary}(\mathrm{x}))}`$ zaksJSP97 , in analogy to eq. (14), the averaged exponents of the resulting spectra are different.

To conclude, the one-dimensional Sierpinski automaton generates $`1/f^\alpha `$ spectra in the number of active states, and can therefore be considered as one of the simplest models generating $`1/f^\alpha `$ spectra. While the Sierpinski automaton is rather a caricature, the approach of studying the sum signal, or total (in-)activity, and its spectrum, can be transferred to more realistic models and compared directly with experiments. Although exact Sierpinski patterns with long-range correlations remain to be experimentally challenging, we conjecture that $`1/f^\alpha `$ spectra in a suitable sum signal can be identified in every experimental setup exhibiting Sierpinski patterns.
Acknowledgements: We would like to thank P.J. Plath, O. Rudzick, and S. Bornholdt for fruitful discussions.


text:
# Ballistic Spin Injection from Fe into ZnSe and GaAs with a (001), (111), and (110) orientation

## I Introduction

The idea of exploiting both the charge and the spin of an electron in semiconductor (SC) devices has lead to the growing field of spintronics. WA01  Several spintronic devices are already proposed, DD90 ; GS99  but the full potential of spintronics has still to be discovered. At the moment the essential ingredients for spintronics, e.g., the injection and detection of a spin polarized current, the spin transport in SC etc., are only partially achieved in experiments. In this paper we want to investigate one of these main challenges: How to create a spin polarized current in a SC? There are several methods that obtain a high spin polarization, but most are useful for basic studies rather than applications: illumination of the SC with circularly polarized light, HO98 ; KA98  injection of electrons from a ferromagnetic scanning tunneling microscope (STM) tip AR92 ; LB01  or from a paramagnetic semiconductor polarized in an external magnetic field.FK99 ; JP00  Better applicable for devices is the injection from a diluted magnetic semiconductor,OY99 ; Mattana03  but these have the drawback of low Curie temperatures $`T_C`$ up to now. More promising is the injection from a common ferromagnetic metal, like Fe, into a SC due to the high $`T_C`$. In this case it has been shown SF00  that in the diffusive limit the spin polarization of the injected current is vanishingly small, FH00 ; GS99 ; HB99  if the contact between the FM and the SC is ohmic. This ‘fundamental obstacle’ is traced back to the conductivity mismatch between both materials. It can be overcome by either a ferromagnet with nearly 100% spin polarization at the Fermi energy, e.g., a half metal, or by inserting a spin polarized interface resistance at the FM/SC interface, Ra00 ; FJ01  e.g., an extrinsic tunneling barrier or an intrinsic Schottky barrier. By means of a Schottky barrier ZR01 ; HJ02 ; HE03  or of an Al2O3 barrier MS02  spin injection could be achieved in experiments.

Here we investigate the spin injection from a ferromagnet (FM) into a SC in the ballistic regime, i.e., the limit where inelastic and incoherent scattering events are negligible. In this case there is no scattering in the perfect bulk regions and the whole resistance results only from the interface region so that there is no conductivity mismatch. Therefore other effects emerge that are not present in the diffusive limit: Kirczenow Ki01  has shown that some FM/SC interfaces can act as ideal spin filter. This is the case, if for a direct-gap SC the Fermi surface of the FM projected onto the two-dimensional Brillouin zone (2DBZ) has a hole at the $`\overline{\mathrm{\Gamma }}`$ point for one spin direction. But the material combinations Fe/ZnSe and Fe/GaAs investigated here, do not fulfill this condition. It has been pointed out that the symmetry of the FM wave functions is very important for the spin injection process. WM02 ; MW02 ; ZX03  If there is a second FM lead to detect the spin polarization of the current, additional effects like the forming of quantum well states, Fabry-Perot like resonances and extremely high magnetoresistance ratios have been reported in a recent theoretical ab-initio study. MW02

There are two types of calculations for the ballistic spin injection already published: by means of analytical models using plane waves Gr01 ; HM01 ; HS01  and by ab-initio methods. WM02 ; ZX03 ; MW02 ; Wu03  While the model calculations obtain only a few percent of spin polarization, the inclusion of the full band structure of the FM and the interface region in ab-initio methods can result in spin polarizations up to 99%. The origin of this high polarization are the different symmetries of the Fe $`d`$ states at the Fermi energy for the majority and the minority spin.

Here we extend our previous work WM02 ; MW02  and show more details of the spin injection process. In Sec. II the details of the calculations and the investigated Fe/SC systems are given. In Sec. III the ballistic spin injection of hot electrons for the Fe/SC(001), (111), and (110) interfaces is investigated. It emerges that the (111) and (110) orientations do not show the large symmetry enforced spin polarization as the (001) orientation does. In Sec. IV we report on injection of thermal electrons with and without a Schottky barrier, addressing also the effect of resonant interface states. Sec. V gives an approximation to the effect of a higher bias voltage (out of the linear-response regime). The paper closes with a summary of the results in Sec. VI.

## II Calculation details

The investigated heterostructures consist of a half-space of ideal bulk Fe and one of the ideal bulk SC (either ZnSe or GaAs). Their properties are described by the surface Green’s function determined by the decimation technique. TD97  Between the half-spaces there are several monolayers (MLs) of the interface region, where the atomic potentials are allowed to deviate from their bulk values. This interface region consists of four MLs of Fe and two MLs of SC for the (001) orientation and appropriate numbers of MLs for the (111) and (110) orientations. The corresponding potentials are obtained self-consistently from a Fe/SC/Fe junction geometry by the screened Korringa-Kohn-Rostoker (KKR) Green function method.WZ97 ; Papanikolaou02  It is assumed that the SC lattice is matched to the bulk Fe, resulting in an fcc SC lattice constant double the lattice constant of bcc Fe $`d_{\mathrm{SC}}=2d_{\mathrm{Fe}}=5.742`$Å. The experimental lattice constant of Fe is used. This means that the SC is stretched compared to the bulk lattice constant by 1.3% and 1.6% for ZnSe and GaAs, respectively, resulting in a slightly smaller energy gap. In order to describe the SC accurately, two empty spheres per unit cell are introduced to account for the open space in the geometry of the zinc-blende structure. A perfect two-dimensional translation symmetry in the whole system is assumed and so the in-plane component $`𝐤_{}`$ of the $`𝐤`$ vector is conserved while crossing the interface. The $`z`$ axis is always assumed to be the growth direction, i.e., standing perpendicular on the interface. Also the conductance $`G`$ is calculated by a two-dimensional Fourier transformation in dependence of $`𝐤_{}`$.

The calculations of the ground state properties are performed within the density-functional theory in the local density approximation (LDA) for the exchange and correlation terms. It is well known that it underestimates the energy gap in SC by around one half. Otherwise it gives very accurate FM and SC bands and is well suited for the considered problem. The potentials are described within the atomic-sphere approximation, and the wave functions are expanded in angular momentum up to a cutoff of $`\mathrm{}_{max}=2`$ for self-consistency and $`\mathrm{}_{max}=3`$ for the conductance calculations. Spin-orbit coupling and spin-flip scattering are not included in the calculations. This could have an effect in case of injection into the SC valence band, where the spin-orbit interaction is stronger, but not in the case of the valence band which has mainly $`s`$ character and negligible spin-orbit coupling.

The conductance $`G`$ is calculated by the Landauer formula. La57  $`G`$ is expressed as a sum of the transmission probabilities through the scattering region (here the Fe/SC interface) over all available conducting channels. These channels are enumerated by the spin orientation $`\sigma =`$ or $``$, the bands $`\nu ,\nu ^{}`$ at the Fermi energy of the left and right lead, respectively, and the $`𝐤_{}`$ vector. The spin-dependent conductance, normalized to the area of the two-dimensional unit cell, reads

$$G^\sigma =\frac{e^2}{h}\frac{1}{A_{2\mathrm{D}\mathrm{B}\mathrm{Z}}}\underset{\nu ,\nu ^{}}{}\underset{2\mathrm{D}\mathrm{B}\mathrm{Z}}{}d^2k_{}T_{\nu ,\nu ^{}}^\sigma (𝐤_{})$$
(1)
with $`A_{2\mathrm{D}\mathrm{B}\mathrm{Z}}`$ the area of the two-dimensional Brillouin zone. Since the $`𝐤_{}`$ vector and the spin $`\sigma `$ of the electron are conserved in our calculation, they are the same for the incoming and the outgoing states. Here the Landauer formula is evaluated by a Green’s function formalism introduced by Baranger and Stone. BS89  Details are presented elsewhere. MPxx  The spin polarization $`P`$ of the current is defined by the spin dependent conductances in both spin bands

$$P=\frac{G^{}G^{}}{G^{}+G^{}},$$

where $`G^{}`$ and $`G^{}`$ is the conductance of the majority and the minority electrons, respectively.

## III Hot electron injection process

To investigate the effects of the symmetry of the Fe bands, we calculate first the hot injection of electrons, i.e., of Fe electrons well above the Fermi level $`E_F`$, which falls in the gap of the SC. For most of the calculations we restrict the discussion to the $`\overline{\mathrm{\Gamma }}`$ ($`𝐤_{}=0`$) point for simplicity, i.e., we consider only electrons with perpendicular incidence on the interface. This is motivated by the fact that in most spin injection experiments the Fermi energy in the SC is only some tens of meV above the conduction band minimum resulting in a very small Fermi sphere around the $`\mathrm{\Gamma }`$ point.

### III.1 The Fe/SC(001) orientation

First we investigate the spin injection in Fe/SC junctions grown in the (001) orientation. The (001) direction has the highest symmetry and is therefore the best candidate for a symmetry enforced high spin polarization of the injected current,WM02  i.e., almost total reflection of the minority spin electrons due to the symmetry mismatch between the SC and minority spin Fe bands.

For the investigation of the bands that are available in Fe and in the SC at the $`\overline{\mathrm{\Gamma }}`$ point, it is important to note that the in-plane common two-dimensional unit cell is double in area than the one of bulk bcc Fe (containing two inequivalent Fe atoms), and the common 2DBZ has half the area of the one of bcc Fe. Thus additional bands are backfolded into the 2DBZ and the number of bands at $`\overline{\mathrm{\Gamma }}`$ increases. A description of the backfolding is given in the Appendix.

The Fe band structure at the $`\overline{\mathrm{\Gamma }}`$ point along $`k_z`$ ($`\mathrm{\Gamma }`$-H, i.e. the $`\mathrm{\Delta }`$ direction), including the backfolded bands, is shown in Fig. 1 for both spins. This direction has in bulk Fe a $`C_{4v}`$ symmetry, just as the Fe (001) surface.

In the 2DBZ of the SC no bands are backfolded to the $`\overline{\mathrm{\Gamma }}`$ point. The ZnSe and GaAs band structures along $`\mathrm{\Gamma }`$-X (the $`\mathrm{\Delta }`$ direction) are shown in Fig. 2. In the SC this direction has a $`C_{2v}`$ symmetry (a subgroup of $`C_{4v}`$), as the SC (001) surface.

Therefore the Fe/SC(001) interface has also the reduced $`C_{2v}`$ symmetry. As a result, some bands $`d`$ that are mutually orthogonal in bulk Fe are able to couple to each other at the interface.

Now it is important to examine which bands of Fe are allowed to couple to the $`\mathrm{\Delta }_1`$ conduction band of the SC which transforms fully symmetrically under the rotations of $`C_{2v}`$; this can be found by the compatibility relations between the groups $`C_{4v}`$ and $`C_{2v}`$. Note that the symmetry notation of Fe and the SC refers to the different groups, $`C_{4v}`$ for Fe and $`C_{2v}`$ for the SC. At $`E_F`$ and above the available Fe bands along the $`\mathrm{\Delta }`$ direction (black lines in Fig. 1) and their symmetry and coupling properties are given in Table 1. Practically, the only ones that are able to couple well are the $`\mathrm{\Delta }_1`$ bands containing the extended $`s`$ and $`p_z`$ orbitals, as well as the $`d_{z^2}`$ orbitals pointing into the SC, while the $`\mathrm{\Delta }_2^{}`$ bands couple only poorly because they consist of more localized $`d_{xy}`$ orbitals pointing in-plane. At $`E_F`$, and up to 1.2 eV above, the former bands exist only for majority spin, allowing transmission, while for minority spin the latter bands co-exist with the $`\mathrm{\Delta }_2`$ and $`\mathrm{\Delta }_5`$ bands.

In Fig. 3 the calculated conductance for the hot injection process for Fe/ZnSe(001) and in Fig. 4 for Fe/GaAs(001) are shown. In (a) the Zn and Ga terminated and in (c) the Se and As terminated interfaces are shown. The intermixed interfaces (b) will be discussed later.

As it was already discussed above, the coupling of the minority $`\mathrm{\Delta }_2^{}`$ band is much weaker than the coupling of the majority $`\mathrm{\Delta }_1`$ band to the $`\mathrm{\Delta }_1`$ conduction band of the SC. Thus the $`\mathrm{\Delta }_2^{}`$ states are nearly totally reflected at the interface. This results for all three terminations in a very high symmetry enforced spin polarization. Around 1.3 eV above the Fermi energy (at 11.3 eV) also a $`\mathrm{\Delta }_1`$ band is available in the minority channel. At this energy the conductance in the minority channel rises drastically showing that the high spin polarization at lower energies originates precisely from the absence of the $`\mathrm{\Delta }_1`$ Fe band in the minority channel. Also the flat backfolded bands of D1 and D2 symmetries give a small contribution to the conductance in the small energy window of around 0.3eV width, centered around $`11.4`$eV. At these energies in the ZnSe half-space there is only one conduction band available at the $`\overline{\mathrm{\Gamma }}`$ point, limiting the maximum conductance to 1 $`e^2/h`$. In the case of the GaAs half-space there are three conduction bands and so the conductance could rise up to a maximal value of 3 $`e^2/h`$.

Finally the conductance can be seen to decrease for small conduction band energies and even vanishes at the conduction band minimum of the SC. This can be explained qualitatively by the vanishing group velocity, so that no current can be transported away from the interface.

### III.2 The Fe/SC(111) interface

In this section the spin injection for the Fe/SC junctions grown in the (111) orientation is investigated. In the Fe half-space the  direction has a hexagonal symmetry with one atom per unit cell. But in the SC half-space this direction has only a three-fold rotational axis. In each ML there is either a cation or an anion (or one of the two types of empty spheres). Also there are a two kinds of possible geometric terminations: one where the terminated SC atom is singly and one where it is triply coordinated to the Fe atoms.

In Fig. 5 the band structure of Fe(111) at the $`\overline{\mathrm{\Gamma }}`$ point in the  direction is shown.

In the SC(111) 2DBZ no bands are backfolded to the $`\overline{\mathrm{\Gamma }}`$ point and the only available band corresponds to the bulk band along the $`\mathrm{\Gamma }`$-L high symmetry line.

The band structure plot for Fe(111) shows important differences compared to the (001) orientation: At the Fermi energy in both spin channels there are bands with the same $`\mathrm{\Lambda }_1`$ symmetry and there cannot be any symmetry enforced spin polarization. All eventually obtained differences in the conductances for both spin directions are due to the different coupling of the Fe bands to the SC. This coupling is believed to be sensitive to the exact interface properties like the termination, lattice relaxations etc. The calculated results for the hot injection process in Fe/ZnSe(111) with a Zn terminated interface are shown in Fig. 6.

Also for this orientation the conductance is shown at the $`\overline{\mathrm{\Gamma }}`$ point as before.

In both the singly and the triply-coordinated Zn terminations the obtained spin polarization is much smaller than in the (001) interface, demonstrating the lack of a symmetry enforced spin filtering. Above 12.0eV there is no ZnSe band exactly at the $`\overline{\mathrm{\Gamma }}`$ point in the  direction, and so no propagating states are available.

### III.3 The Fe/SC(110) interface

This section deals with Fe/SC(110) oriented junction. This is the one with the lowest symmetry considered in this work. From this point of view a high symmetry-enforced spin polarization is most improbable. Also the surface unit cell in the SC half-space contains all four kinds of atoms (cation, anion and both types of vacancies) and is larger than the (001) and the (111) unit cells. By matching the two dimensional translational symmetry, also the Fe surface unit cell has to contain four Fe atoms. This results in a very small 2DBZ for this orientation leading even in the SC half-space to backfolded bands at the $`\overline{\mathrm{\Gamma }}`$ point.

The band structure of Fe (001) at the $`\overline{\mathrm{\Gamma }}`$ point with the different backfolded bands are shown in Fig. 7.

In the SC(110) half-space there is also a backfolded band at the $`\overline{\mathrm{\Gamma }}`$ point (created by applying $`𝐛_1`$ in Fig. 18). Because of this the maximum conductance per spin direction is $`2e^2/h`$.

The variety of backfolded bands and the reduced symmetry of the  direction do not allow for a simple symmetry-based discussion. Basically, for the (110) oriented heterojunctions there is no symmetry enforced spin polarization. This is also shown in Fig. 8, where the conductances at the $`\overline{\mathrm{\Gamma }}`$ point for the hot injection process are presented.

From the above we conclude that the best candidate for a symmetry-enforced high current spin polarization is the (001) interface. Thus, in what follows, we restrict our study to this case.

### III.4 Hot injection process $`𝐤_{}`$ resolved

In the hot injection process the restriction to the $`\overline{\mathrm{\Gamma }}`$ point was motivated by the fact that in most applications the Fermi level is positioned only slightly in the conduction band resulting in a very small Fermi sphere around the $`\overline{\mathrm{\Gamma }}`$ point. More correctly the conductance should be integrated over the whole 2DBZ, because also states away from the $`\overline{\mathrm{\Gamma }}`$ point are populated for higher injection energies. Therefore in Fig. 9 the $`𝐤_{}`$ resolved conductances are shown for the majority and minority electrons for the Zn-terminated Fe/ZnSe (001) interface.

The energy of the injected electrons is 500 meV above the conduction band minimum. In this case the polarization at the $`\overline{\mathrm{\Gamma }}`$ point is 99.7%. If the conductance is integrated over the 2DBZ the polarization is reduced to 39.5%. The reason is a higher conductance in the minority band for states away from the $`\overline{\mathrm{\Gamma }}`$ point. The symmetry arguments discussed above are valid only at this high symmetry point. For states $`𝐤_{}0`$ also other Fe minority states are allowed to couple to the SC $`\mathrm{\Delta }_1`$ conduction band and so the conductance rises. We conclude that the high symmetry-enforced spin polarization can be realized only for energies up to a few tens of meV above the conduction band edge.

In the conductance plot (Fig. 9) for the minority electrons the reduced $`C_{2v}`$ symmetry of the Fe/SC(001) interface can be directly seen. The majority conductance behaves like free electrons traveling across a potential step and is practically unaffected by the reduced symmetry. The circularly symmetric form of $`G^{}(𝐤_{})`$ and the twofold-symmetric form of $`G^{}(𝐤_{})`$ can be understood in terms of the form of the Bloch functions for small deviations from $`𝐤_{}=0`$ in the same manner as described in Ref. MW02,  for the Fe/SC/Fe (001) junctions.

### III.5 The intermixed interface

By means of an ab-initio study EL02  it has been found that the interface of Fe/GaAs(001) with a Ga terminated interface is energetically more stable, if Fe atoms diffuses into the vacancy sites of the first Ga ML. In case of an As terminated interface the abrupt interface is found to be more stable than the intermixed one. The former case is accompanied by rather large lattice distortions which are not taken into account in our calculations. Although we do not know of analogous studies for Fe/ZnSe(001) interfaces we have also performed the intermixed interface in case of a Zn termination. The results are shown in Fig. 3(b) and 4(b) for the hot injection process with an intermixed Zn and Ga interface, respectively. The still high spin injection polarization can be explained by the fact that the intermixed Fe atoms do not change or reduce the $`C_{2v}`$ symmetry of the interface and so the symmetry arguments are still valid. This is a different effect than interface roughness that reduces the spin polarization drastically. ZX03 ; SH02  Solely the coupling of the Fe states to the SC are changed. In both cases the coupling in the minority band is enhanced in comparison to the atomically abrupt interfaces shown in Fig. 3(a) and 4(a).

## IV Thermal injection process

In this section the injection of thermal electrons is investigated. This is achieved by lowering the potentials in the SC half-space by a rigid shift, so that the Fermi energy falls slightly above the conduction band minimum (here around 10 meV). All other potentials are kept in their ground state position, in particular also the potentials of the two SC monolayers in the interface region. This shifting simulates the effect of n-doping or an applied gate voltage in a field-effect-transistor device. The electrons are injected at the $`E_F`$ of the Fe half-space directly into the conduction band. Since the Fermi energy is only some tens of meV above the conduction band minimum and since effective masses of ZnSe and GaAs are small, the resulting Fermi wave vector in the SC half-space is very small (around one hundredth of the distance to the boundary of the 2DBZ). So we present here the conductance only for the $`\overline{\mathrm{\Gamma }}`$ point. The evaluated polarization changes only slightly, if the conductances are integrated over the whole 2DBZ.

### IV.1 Thermal injection without Schottky barrier

The results for Fe/ZnSe and Fe/GaAs(001) are discussed in Ref. WM02,  and show a very high spin polarization of more than 97% for ZnSe and practically 100% for GaAs (001). Similar results are found for tunneling in the Fe/SC/Fe (001) junctions.MW02  Here we add the thermal injection with an intermixed Zn- and Ga-terminated interface (Fig. 10).

In the case of an intermixed Zn interface nearly the same conductances are obtained as for the abrupt one (Fig. 3 in Ref. WM02, ). But a drastic increase of the minority conductance can be seen for the intermixed Ga interface, resulting in a much smaller spin polarization. The reason for this increase is an interface resonance localized at the intermixed Fe atoms. This can be seen in the density-of-states (DOS) for this heterojunction shown in Fig. 11 for the intermixed Zn and Ga interface. The arrows indicate the position of the Fe resonance. In the case of the intermixed Zn interface this resonance lies slightly above the $`E_F`$ and give rise to the small maximum in the minority conductance at $`E_FE_C=15`$meV seen in Fig. 10.

In the DOS for the intermixed Ga interface this resonance lies nearly on the Fermi energy and is much narrower, so that it has a strong influence on the minority conductance. Since the energy position of this resonance is believed to depend strongly on the interface properties, e.g., lattice relaxations, further investigations are needed for the exact energy position. If the injection process is evaluated at only 0.3eV higher energies, the influence of this resonance is negligible and the normal high spin polarization is restored. However, for a higher bias one must integrate over an energy interval (see below) and thus include resonant interface states that lie close to, but not at, $`E_F`$.

### IV.2 Thermal injection through a Schottky barrier

In theoretical studies for strongly diffusive transportRa00 ; FJ01  it is shown that the negligible spin polarization due to the conductivity mismatch between the ferromagnetic metal and the semiconductorSF00  can be overcome by a spin dependent tunneling barrier at the interface. Such a barrier could be the Scottky barrier created at the interface. Although the approaches in Refs. Ra00, ; FJ01, ; SF00,  assume diffusive transport, tunneling barriers can be described in the ballistic regime.

In this section the thermal injection process is extended to the case where a Schottky barrier at the interface is included. The potentials for this junctions are the same as in the thermal injection process, but with an inserted tunneling barrier at the interface. The barrier is modeled by a rigid shift of each SC atom potential, constant within each monolayer. It starts at the third SC layer from the interface, and increases linearly in magnitude with distance from the interface. In this way the SC bands are gradually shifted downward. Finally, at the end of the barrier $`E_F`$ lies 10 meV above the conduction band edge (and the shift is not changed from then on), while at the interface $`E_F`$ lies in the middle of the gap. The barrier thickness is then varied between 0 and 140 MLs.

As shown in Ref. WM02,  a tunneling barrier at the Fe/SC (001) interface gives also a high spin polarization. Only in case of a resonant interface state near $`E_F`$, especially in Fe/ZnSe(001) with a Zn terminated interface, the polarization is reduced. Here we will discuss the effect in more detail. In Fig. 12 the conductance and the polarization in dependence of the barrier thickness is shown for Fe/ZnSe and Fe/GaAs(001), respectively.

Due to the small Fermi wave vector, the conductance is again analyzed only at the $`\overline{\mathrm{\Gamma }}`$ point. The polarizations agree with the ones correctly integrated over the 2DBZ within 5%.

Except for the Se terminated Fe/ZnSe(001) interface, all terminations show the influence of a resonant interface state in the minority band, resulting in an higher conductance for the minority spins with thicker barriers. This leads to a reduction of the spin polarization. For Fe/GaAs(001) this effect is also visible in the conductance but is negligible for the spin polarization due to the much larger difference of the conductance for both spins than in Fe/ZnSe. To discuss this effect in more detail, in Fig. 13 the DOS at the $`\overline{\mathrm{\Gamma }}`$ point is shown for Fe/ZnSe(001). The DOS for Fe/GaAs(001) is qualitatively comparable for this purpose.

In the DOS an interface state is visible near the Fermi energy in the minority spin channel. This state has $`\mathrm{\Delta }_1`$ symmetry and lies in the energy region of the $`\mathrm{\Delta }_1`$ hybridization gap of the Fe minority band structure. It can penetrate well a Schottky barrier of moderate thickness since the evanescent wave with the longest decay length is of the same $`\mathrm{\Delta }_1`$ symmetry.MP00  In the Fe half-space this interface state is normally also evanescent since there are no propagating states it can couple to. But due to the reduced symmetry of the interface the $`\mathrm{\Delta }_1`$ interface state couples weakly to the Fe $`\mathrm{\Delta }_2^{}`$ band and becomes resonant. In the case of the Se terminated Fe/ZnSe(001) interface the interface state lies further away from $`E_F`$ and practically does not contribute to the conductance.

The potentials used for the DOS in Fig. 13 are the ground state potentials, i.e., $`E_F`$ is in the middle of the gap in the SC half-space and no Schottky barrier is inserted. If now a Schottky barrier is introduced into the junction, the system and the resonant interface state are slightly distorted, because the $`\mathrm{\Delta }_1`$ interface state interacts with the $`\mathrm{\Delta }_1`$ conduction band of the semiconductor. This results in a slight downshift in energy of the resonant interface state away from $`E_F`$ for smaller Schottky barriers. The reduction of the spin polarization with thicker Schottky barriers can thus be explained by a shift of the resonant interface state towards $`E_F`$.

This effect can also be seen by changing the pinning position of the Fermi energy relative in the SC gap at the interface (the potentials in the interface region are kept fixed). In Fig. 14 the conductance for different pinning positions in the gap is shown for Zn-terminated Fe/ZnSe(001) with a 80 ML thick Schottky barrier. The minority-conductance peak at $`0.2`$ eV reflects the influence of the interface state.

Furthermore one can estimate the interface resistance under the assumption of diffusive transport in the bulk, using Schep’s formula.SH97 ; SP00  Although this has been derived for metallic multilayers, we have applied it to this case. Our results are summarized in Table 2. More details can be found in Ref. Wunnicke03, .

## V Higher bias voltage

Up to now it was always assumed that the applied bias voltage is so small that only states at $`E_F`$ carry current. But experimentally especially for thicker barriers, a bias voltage in the order of 1V or less is used. HJ02  This relatively high bias in experiments with an optical detection of the spin polarization is needed to gain a reasonable signal-to-noise ratio of the emitted light. It is believed that the bias can be reduced in experiments with thinner tunneling barriers or with an electrical detection of the polarization, since then a much smaller current is needed. If a non-zero bias is applied, one must integrate the transmission over the appropriate energy range, and the Landauer formula readsCV97

$`G={\displaystyle \frac{1}{eU}}{\displaystyle \underset{0}{\overset{eU}{}}}𝑑E`$ (2)
$`{\displaystyle \frac{e^2}{h}}{\displaystyle \frac{1}{A_{2\mathrm{D}\mathrm{B}\mathrm{Z}}}}{\displaystyle \underset{\sigma ,\nu ,\nu ^{}}{}}{\displaystyle \underset{2\mathrm{D}\mathrm{B}\mathrm{Z}}{}}d^2k_{}T_{\nu ,\nu ^{}}^\sigma (𝐤_{},E_F+E)`$
with $`T_{\nu ,\nu ^{}}^\sigma (𝐤_{},E_F+E)`$ being the transmission probability at the energy $`E_F+E`$. This equation gives the standard Landauer formula in the case of small bias voltages $`U`$ in a first approximation, under the assumption that the electronic structure is unaffected by the applied voltage.

Here we show that also for a non-zero bias and a tunneling barrier of moderate thickness the $`\overline{\mathrm{\Gamma }}`$ point is of main importance and a high spin polarization is obtained. It has been shown MP00  that in tunneling the smallest damping factor $`\kappa `$ is located at the $`\overline{\mathrm{\Gamma }}`$ point in the energy gap for ZnSe and GaAs (and for other direct-gap semiconductors), so that this point plays the major role for thick tunneling barriers. Also at this point the symmetry enforced spin polarization is the highest, as shown in the $`𝐤_{}`$-resolved hot injection process. In the calculations the Fermi energy in the bulk SC is assumed to be 10 meV above the conduction band minimum. For the Zn termination an abrupt interface is taken for simplicity. Due to the high bias voltage a $`𝐤_{}`$ integration over the whole 2DBZ is performed, because also states away from the $`\overline{\mathrm{\Gamma }}`$ point in 2DBZ are populated and carry the current.

The position of the energy bands and of the Fermi energy in the Fe and SC half-spaces are sketched in Fig. 15.

Already from this simple sketch it is qualitatively clear that the conductance at the upper energy range (at the Fermi energy of the Fe half-space) is the most important one, because there the effective tunneling barrier is the smallest. In Table 3 the spin polarizations are listed for a Fe/ZnSe(001) junction with an applied voltage of 1V and a tunneling barrier of 8 ML and 80 ML thickness.

In case of a Zn termination it is important that slightly below the Fermi energy of the Fe half-space the resonant interface state contributes largely to the minority conductance. This gives a relatively low and a more or less constant spin polarization with thicker tunneling barriers. In the Se terminated interface the resonant interface state lies lower in energy, where the effective tunneling barrier is thicker than for the Zn terminated interface. Thus the resonant interface state has a much smaller influence on the minority conductance. With larger barrier thickness the high spin polarization is more or less restored because only states at $`\overline{\mathrm{\Gamma }}`$, having the smallest decay parameter are important and the symmetry arguments apply again.

## VI Summary and conclusions

We have reported first-principles calculations on spin injection from Fe into GaAs and ZnSe through the (001), (111), and (110) interfaces. The electronic transport has been assumed to be in the ballistic regime, and the interfaces have been assumed to have two-dimensional periodicity, so that $`𝐤_{}`$ is conserved during scattering at the interface. Then, in a Landauer-Büttiker approach, the conductance is determined by the transmission probability summed for all $`𝐤_{}`$. Under these assumptions, we have considered hot and thermal injection — the latter also in the presence of a Schottky barrier — and approached the higher bias regime. We have reached the following conclusions.

(i) The spin polarization of the current is highest when most incoming Fe Bloch states of one spin direction are totally reflected due to symmetry mismatch with the SC conduction band states — we call this effect “symmetry-enforced spin polarization”. This is the case for the (001) interface, where minority-spin electrons for $`𝐤_{}=0`$ practically cannot be transmitted. The lower symmetry of the (111) and (110) interface does not lead to such a selection rule. Even for the (001) interface, in the case of hot injection, above an energy threshold where the Fe minority $`\mathrm{\Delta }_1`$ band starts, the minority-spin current increases rapidly even at $`𝐤_{}=0`$, being then symmetry-allowed.

(ii) For the symmetry-enforced spin polarization to be realised, the injection must take place close to the conduction band edge, so that only states near $`𝐤_{}=0`$ are relevant; these have a well-defined and suitable symmetry properties. Also, the interface should be as ordered as possible, otherwise $`𝐤_{}`$ is not conserved and incoming minority-spin Fe states from all $`𝐤_{}`$ can scatter into the SC conduction band. Then the spin polarization will decrease.

(iii) The case of thermal injection (i.e. exactly at $`E_F`$), has been considered for the (001) interface with and without a Schottky barrier. The symmetry concepts still hold, since the least-decaying complex band of the tunneling barrier has the same symmetry properties as the conduction band (both are of $`\mathrm{\Delta }_1`$ character). We have seen that resonant interface states existing in the vicinity $`E_F`$ for minority spin can be important for the current spin polarization. Even in the case of a Schottky barrier they can provide a resonant tunneling channel for the minority-spin electrons decreasing the spin injection effect.

(iv) In the case of a higher bias, a wider energy range must be considered, so that states with $`𝐤_{}0`$ as well as resonant states close enough to $`E_F`$ become relevant. Then the spin polarization decreases, but spin injection is still achieved.

It can be argued that a well ordered interface, necessary for our symmetry selection rule, is difficult to realize experimentally. However, recent successful attempts in increasing the spin injection efficiency have been reported, accompanied by an improvement of the quality of the interface.HE03  Therefore, we are optimistic that our work will motivate further research in this direction.

###### Acknowledgements.

The authors thank Arne Brataas for helpful comments. This work was supported by the RT Network Computational Magnetoelectronics (Contract RTN1-1999-00145) of the European Commission.

## Appendix A Geometry of the reduced 2DBZ in (001), (111), and (110)

First we investigate the backfolded bands for the Fe(001) half-space. In Fig. 16 the 2DBZ and the cut perpendicular to the  direction of the bulk Brillouin zone are shown.

At the $`\overline{\mathrm{\Gamma }}`$ point two bands are available: the not backfolded band corresponding to the bulk band along $`\mathrm{\Gamma }`$-H, i.e., the $`\mathrm{\Delta }`$ direction. This direction has a $`C_{4v}`$ symmetry. Also one band is backfolded by applying once the reciprocal surface lattice vectors $`𝐛_\mathrm{𝟏}`$ or $`𝐛_\mathrm{𝟐}`$, corresponding to the bulk band along N-P-N, i.e., the D direction.

Next we discuss the (111) case. In Fig. 17 the (111) cut through the bulk Brillouin zone and the 2DBZ are shown. As for the (001) orientation additional bands are obtained by backfolding to the $`\overline{\mathrm{\Gamma }}`$ point due to the large lattice constant in the Fe half-space.

The not backfolded band at the $`\overline{\mathrm{\Gamma }}`$ point corresponds to the bulk band along the $`\mathrm{\Gamma }`$-P-H high symmetry line, because in the 2DBZ the $`k_z`$ vector extends from $`\sqrt{3}`$ to $`\sqrt{3}`$ in units of $`2\pi /d_{Fe}`$. By applying $`𝐛_1`$ or $`𝐛_2`$, together one additional band is backfolded to the $`\overline{\mathrm{\Gamma }}`$ point. It corresponds to the bands of the same $`\mathrm{\Gamma }`$-P-H high symmetry line as the not backfolded band.

The 2DBZs of Fe and SC(110) and the corresponding (110) cut through the bulk Brillouin zones are shown in Fig. 18.

In Fe(110) there are four backfolded bands by applying the reciprocal surface lattice vectors $`𝐛_1`$ and $`𝐛_2`$. The not backfolded band corresponds to the $`\mathrm{\Gamma }`$-N bulk band ($`\mathrm{\Sigma }`$ direction). By applying the $`𝐛_1`$ reciprocal lattice vector once, one ends up with a backfolded band that is not along a high symmetry line. More precisely, it is along the shortest line that connects two P points in the bulk Brillouin zone. Written in components of the reciprocal bulk lattice vectors, this backfolded band corresponds to the bulk band along

$$\left(\begin{array}{c}1/4\\ 1/4\\ 1/2\end{array}\right)\frac{2\pi }{d_{\mathrm{Fe}}}\left(\begin{array}{c}1/4\\ 1/4\\ 1/2\end{array}\right)\frac{2\pi }{d_{\mathrm{Fe}}}.$$

If both reciprocal lattice vectors, $`𝐛_1`$ and $`𝐛_2`$, are applied, the backfolded band is parallel to the one, if only $`𝐛_1`$ is applied, but shifted in the (110) direction. Written in components of the reciprocal bulk lattice vectors, it corresponds to

$$\left(\begin{array}{c}1/4\\ 1/4\\ 1/2\end{array}\right)\frac{2\pi }{d_{\mathrm{Fe}}}\left(\begin{array}{c}3/4\\ 3/4\\ 1/2\end{array}\right)\frac{2\pi }{d_{\mathrm{Fe}}}.$$

The next backfolded band is obtained by applying $`𝐛_2`$ once that corresponds to the bulk band between the high symmetry points N-H. In summary due to the large two-dimensional surface unit cell and small 2DBZ there are four bands available at the $`\overline{\mathrm{\Gamma }}`$ point for the (110) orientation in the Fe half-space including two bands corresponding *not* to a bulk band along a high-symmetry line. For these bands the only symmetry operation is the identity one, so that they are allowed to couple to any SC bands due to symmetry.

In the SC (110) SBZ the not backfolded band corresponds to the bulk high symmetry line $`\mathrm{\Gamma }`$-$`K`$ ($`\mathrm{\Sigma }`$-direction). By applying $`𝐛_1`$ one additional band is backfolded to the $`\overline{\mathrm{\Gamma }}`$ point, exactly on the not backfolded band. So the maximum conductance can be 2 $`e^2/h`$ per spin channel.


text:
# Electronic structure and magnetism of Mn doped GaN

## I Introduction

Diluted magnetic semiconductors (DMS) are considered to be potential candidates for present and future technological applications in semiconductor spintronics bookdms . During the last decade, there have been numerous experimental and theoretical studies of II-VI, III-V and IV-VI DMS. Among the III-V DMS, Mn doped GaAs system has been studied rigorously for the last few years. This system shows a Curie temperature ($`T_C`$) of 110 K for a Mn doping concentration of 10 $`\%`$. Despite several attempts, $`T_C`$ couldn’t be raised beyond 175 K. Recently, there have been reports of some room temperature DMS. They include Mn doped GaP mngap , Mn doped chalcopyrite CdGeP2 chalco , Mn doped GaN sasaki ; reed  etc. The origin of ferromagnetism in these compounds is still under debate origin .

Dietl et al.dietl  predicted theoretically a high Curie temperature ($``$ 400 K) for Mn doped GaN (5 $`\%`$ Mn). Their theory was based on a mean field model of hole mediated ferromagnetism. As the Curie temperature for a Mn doped GaAs system is comparatively lower, the theoretical prediction for a higher Tc drew much attention. Sasaki et al. sasaki  grew wurtzite Mn doped GaN films by the molecular beam epitaxy method. Magnetic measurements showed a very high Curie temperature of about 940 K. They ruled out the possibility of phase segregation of some ferromagnetic compound e.g. MnGa and Mn4N which also have high Curie temperatures. Room temperature ferromagnetism in Mn doped GaN was also observed by Reed et al. reed . Deep level optical spectroscopy measurements korotkov  show that Mn forms a deep acceptor level at 1.42 eV above the valence band maximum for GaN doped with small concentration of Mn. It is to be noted that Mn forms an acceptor level at 0.11 eV above the valence band maximum in the case of Mn doped GaAs. So, the overlap of the Mn d-states with the valence band is rather strong in Mn doped GaAs compared to Mn doped GaN.

Recently, there have been a few first principles electronic structure calculations of Mn doped GaN systems. Fong et al.fong  performed electronic structure calculations of Fe and Mn doped GaN using the tight-binding linearized muffin tin orbital (TB-LMTO) method. Sato and Katayama-Yoshida sato  performed KKR-CPA (Korringa-Kohn-Rostoker-Coherent Potential Approximation) calculations to study the relative stabilities of ferromagnetic and spin glass phases. They showed that for a low concentration of Mn, ferromagnetism is favored whereas for the high concentration, the spin glass phase is stable. The disordered local moment model was assumed to describe the spin glass phase. They explained the origin of ferromagnetism in these systems by a competition between double exchange and superexchange interactions. Kulatov et al. kulatov  studied electronic, magnetic and optical properties of zinc-blende (Mn,Ga)N for different concentrations of Mn by the TB-LMTO method in a supercell approach. Anomalous exchange interactions in III-V DMS were found from calculations by Schilfgaarde and Mryasov schilf . They predicted aggregation of magnetic nanoclusters inside the III-V host. Kronik et al. kronik  considered (Mn,Ga)N in wurtzite structure and performed electronic structure calculations using a plane wave pseudopotential method. In a recent preprint alessio , the different origins of ferromagnetism in (Mn,Ga)As and (Mn,Ga)N systems were discussed. The authors pointed out from self-interaction corrected (SIC) pseudopotential calculations that (Mn,Ga)N is characterized by localized Mn 3$`d`$ states with a strong self-interaction. In (Mn,Ga)As, $`d`$ states are weakly correlated and are rather delocalized being strongly hybridized with As $`p`$ states. In this communication, we attempt to understand the electronic structure and magnetic interactions in Mn doped GaN and GaAs systems. The motivation of this paper is twofold : (a) to investigate electronic structure and magnetism of Mn doped GaN system in detail and (b) to have a comparison with Mn doped GaAs system. The paper is organized as follows : In the next section, we describe the computational details. In the results section, a subsection describes the electronic structure and magnetism of the (Mn,Ga)N with varying Mn concentration in the local spin density approximation (LSDA). Then we present results from LSDA+U calculations. Finally, we show a comparison of exchange interactions in (Mn,Ga)N and (Mn,Ga)As systems.

## II Computational details

GaN can be grown both in the zinc-blende and wurtzite structures. But, usually, Mn is doped in a wurtzite GaN host sasaki . In our calculations, we have considered the wurtzite structure. Also, for a comparison, we show calculations for the zinc-blende structure. Experimental lattice parameters such as a=3.189 Å and c=5.185 Å with a c/a ratio of 1.626 were taken for the calculations in the wurtzite structure. Results from the atomic relaxations revealed that the nearest neighbor bond-lengths between Mn and N change only by 3 $`\%`$ compared to that of the bulk GaN. This is in agreement with the results of Kronik et al. kronik . In general, Mn doping in substitutional site results in a small relaxation msm  of the nearest neighbor anions around it.

Calculations have been performed by an ab-initio plane wave code (VASP) vasp . Vanderbilt vander  type ultrasoft pseudopotentials were used for the LSDA calculations. LSDA+U calculations were done in the projector augmented wave (PAW) paw  method as implemented by Bengone et al. bengone . Ga 3d orbitals were included in the basis set of the Ga pseudopotential. A kinetic energy cut-off of 350 eV was used for the plane waves included in the basis set. Ceperley and Alder ca  exchange-correlation functional parameterized by Perdew and Zunger pz  was considered within LSDA. We have also checked that the results obtained within GGA (generalized gradient approximation) gga  are similar. A 8x8x6 k-points grid was used in the Monkhorst Pack scheme mp  for small supercells. For the largest supercell considered, a 2x2x1 grid was used. Local properties such as local density of states and local magnetic moments were calculated by projecting the wave functions onto spherical harmonics eichler . The radii chosen for the projection were 1.31, 1.21 and 0.74 Å for Mn, Ga and N respectively.

We have modeled the system using different supercell sizes to simulate different Mn concentrations. For wurtzite structure, supercells having 4, 8, 16, 32, 72 and 108 atoms were used to model a composition MnxGa1-xN for $`x`$=0.5, 0.25, 0.125, 0.0625, 0.028 and 0.018 respectively. For zinc-blende structure, we used a 64 atom-cell to simulate a Mn concentration of $`x`$=0.03125. It has been found experimentally that Mn occupies the Ga-site expt\_jap . Therefore, in our calculational unit cell, 1 Ga atom was substituted by a Mn atom.

## III Results

### III.1 Calculations within LSDA

In a III-V semiconductor, a cation vacancy creates 3 holes in the valence band leaving anion dangling bonds. When Mn occupies the cation site, it donates 3 electrons to fulfill the bonding. Mn is left with 4 unpaired $`d`$ electrons which give rise to 4 $`\mu _B`$/Mn atom. In a realistic situation, there can be compensating donors e.g. As antisites and interstitial Mn atoms larsprb , present in the system to increase or decrease the magnetic moments. But as we are dealing with ideal uncompensated systems, we always obtain a magnetic moment of 4 $`\mu _B`$/Mn atom for substitutional Mn.

In Fig. 1(a-e), we show the density of states (DOS) for various concentrations($`x`$) of Mn in MnxGa1-xN in the wurtzite structure. In (a), the DOS of the undoped GaN is presented. The large band gap is evident from the figure. The calculated direct band gap at the $`\mathrm{\Gamma }`$ point is 1.9 eV which is underestimated compared to the experimental band gap of $``$ 3.4 eV. This well known underestimation is inherent in the formulation of density functional theory and is well documented in existing literature. From Fig. 1(b) to 1(e), we show the DOSs with increasing $`x`$. The Mn impurity $`d`$ peak in the band gap is away from the valence band for small $`x`$. As the concentration increases, this peak is broadened due to the overlap of Mn $`d`$ wavefunctions and the gap between the impurity peak and valence band edge vanishes. The total DOS/cell gradually regains the shape of the DOS of undoped system as $`x`$ is reduced. However, up to $`x`$=0.25, the Fermi level cuts only the spin up DOSs. As there is no state at Fermi level for the spin down channel, we obtain a half metallic solution giving rise to an integer magnetic moment of 4 $`\mu _B`$/Mn atom. For $`x`$=0.5 shown in Fig. 1(e), the system behaves like a ferromagnetic metal with a reduced magnetic moment. One noticeable difference between Mn doped GaAs and Mn doped GaN systems is the position of Mn impurity band in the gap. In (Mn,Ga)As, Mn impurity band almost merges with the top of the valence band (approximately 0.1 eV above the valence band edge) whereas in (Mn,Ga)N, the Mn impurity band is separated from the valence band by 0.56 eV. The width of this impurity band is 0.94 eV. As the spin polarized Mn impurity band is distinctly in the gap, the valence band is less spin polarized in (Mn,Ga)N compared to (Mn,Ga)As. Spin polarization ($`ϵ_{VBM/CBM}^{}ϵ_{VBM/CBM}^{}`$) of the valence and conduction bands at the $`\mathrm{\Gamma }`$ point is 0.05 eV and -0.42 eV respectively for (Mn0.0625Ga0.9375)N. $`ϵ`$ is the eigenvalue and VBM and CBM represent valence band maximum and conduction band minimum respectively.

In Table I, we present the local charges and magnetic moments of Mn and N atoms for different concentrations of Mn. Charges and magnetic moments of Mn remain almost the same with the concentration variation of Mn. The insensitivity of magnetic moment with concentration is a signature of localized d-states of Mn. The total moment/cell is always 4.0 $`\mu _B`$ which is the signature of a half-metallic solution. The only exception is the case of $`x`$=0.5, where the total moment/cell is 2.77 $`\mu _B`$. The exchange splitting of Mn d-states is less in this case allowing both spin-up and spin-down d-states to cross the Fermi level. The averaged induced moments on nearest neighbor N atoms are also tabulated. In most of the cases, the moments are parallel to Mn moments which is not the case for a Mn doped GaAs system where nearest neighbor As moments are antiferromagnetically coupled to the Mn magnetic moment.

As GaN can be grown in both zinc-blende and wurtzite structures, we have also done calculations for a 6.5 $`\%`$ Mn doped GaN system in zinc-blende structure. In Fig. 2, DOSs for both structures are shown. DOS for a zinc-blende structure is similar to that calculated by Kulatov et al. kulatov . Also the magnetic moment on Mn atom (3.4 $`\mu _B`$) agrees very well. There is no striking difference in the broad features of the DOSs for the two structures. In the zinc-blende structure, the peak at the Fermi level is sharper and the nature of these states are different ($`t_2`$ compared to $`e`$ states for the wurtzite structure). See fig. 5 and related discussions. For both the structures, the Mn impurity peak in the energy gap of the host is separated from the valence band of the host.

### III.2 LSDA+U calculations

It is a matter of debate whether the itinerant band model or the localized atomic model is appropriate for the description of Mn doped semiconductors. Density functional calculations sanvito ; freeman ; msm  based on LSDA or GGA provide the basis of itinerant picture whereas others models are based on a localized atomic picture dietl . A recent photoemission experiment kanski  on Mn doped GaAs system revealed the main Mn d-peak to be situated 3.4 eV below the Fermi level. A previous photoemission experiment okabayashi  reported the peak to be 4.4 eV below the Fermi level. Also, the importance of electron correlation effects in these systems was highlighted. Park et. al. park  argued from LSDA+U calculations that correlation corrections are important to have a better agreement with photoemission spectra. On the other hand, all density functional calculations based on LSDA show a peak around 2.6-2.9 eV below the Fermi level msm ; sanvito . So it can be argued that the completely localized picture or the completely itinerant picture cannot solely describe these systems satisfactorily.

We have done LSDA+U calculations to investigate the effect of electron correlations. Firstly, these calculations do not exist in literature and secondly, it is interesting to compare this simple technique with more rigorous SIC calculations. In Fig. 3(a-d), we show the DOSs obtained from LSDA and LSDA+U calculations for 6.25 $`\%`$ Mn. As the value of U for Mn $`d`$ states is not obtained self-consistently from first principles calculations, we varied U from 4 eV to 7 eV treating it as a parameter. In all cases, the exchange parameter J was considered to be 1.0 eV. Increasing U results in a slow shift of the spin-up impurity band towards the valence band. The small peak around 1.5 eV below the Fermi level ($`E_F`$) is diminished gradually in magnitude whereas the peak around -4.5 eV below $`E_F`$ increases in magnitude. Even for U=10 eV (not shown here), the localized peak remains pinned close to $`E_F`$ and is not merged with the delocalized valence band states. The spin down DOS shifts almost rigidly away from the Fermi level towards higher energy with increasing U. Local charges and magnetic moments of Mn and nearest neighbor N atoms are listed in Table II. Magnetic moment of Mn increases with U due to increase in localization of $`d`$ states. The induced moment of N atoms also increase. In Fig. 3(c), a comparison between LSDA and LSDA+U calculations is shown. The redistribution of weights of Mn $`d`$ peaks with the inclusion of U is visible. In LSDA+U result, the sharp peak is very close to the Fermi level in the impurity band whereas it is in the valence band in the LSDA calculation. This peak is of $`e`$-character and doesn’t take part in the bonding with the neighboring N atoms. These results are in agreement with the more sophisticated SIC calculations reported by Filippetti et al. alessio . They also found a flat band with $`d_{z^2}`$ character at $`E_F`$. The nature of the states close to $`E_F`$ from our calculation is shown later in fig. 5. To our knowledge, the valence band photoemission spectra of (Mn,Ga)N is not available in the literature. So, the extent of validity of LSDA approach cannot be tested. Our LSDA+U findings can be compared with future angle-resolved photoemission experiments on (Mn,Ga)N to verify the existence of the localized peak close to $`E_F`$. A comparative study with Mn doped GaAs is shown in Fig. 3(d). The main broad Mn $`d`$ peak around 2.8 eV below the Fermi level in an LSDA calculation is shifted 4 eV below the Fermi level with a smaller band width. The DOS at the Fermi level is also decreased compared to that of an LSDA calculation. In both LSDA and LSDA+U calculations, the hybridization between Mn $`d`$ and As $`p`$ states are seen. So the holes in the valence band have hybridized $`pd`$ character.

In Fig. 4, we show the integrated magnetization density around various atoms in the unit cell as a function of the radius of integration. It has been calculated as

$$M(R)=_0^R(\rho ^{}(r)\rho ^{}(r))𝑑r$$

where $`\rho ^{}(r)`$ and $`\rho ^{}(r)`$ are the spin-up and spin-down charge densities respectively and $`M`$ is the magnetic moment obtained for a radius $`R`$. For (Mn,Ga)N and (Mn,Ga)As in the zinc-blende structure, the magnetic moment reaches the value 4 $`\mu _B`$ at a distance of 5 Å far from the Mn center. But, in the wurtzite structure of (Mn,Ga)N, this value is reached at a smaller distance. This again shows a more localized character of Mn $`d`$ states in wurtzite GaN. The integration around nearest neighbor anions reveal that N has a positive contribution in (Mn,Ga)N whereas, As has a negative contribution.

In Fig. 5, we show the character of the states within an energy interval close to the Fermi level. The states close to $`E_F`$ in the electronic structure are important in characterizing the origin of ferromagnetism. The projection of the wavefunction onto spherical harmonics around each atom has been calculated as mentioned in ref. eichler . In the figure, $`\chi `$ is defined as : $`\chi _{lm}^{N,\sigma }=_{nk}|Y_{lm}^{N,\sigma }|\varphi _{nk}^\sigma |^2`$, where $`\sigma `$ and $`N`$ are the spin and atom indices respectively and $`n`$ is the index for the bands within the specified energy interval. It is clear from the figure that for (Mn,Ga)N in the wurtzite structure (left panel), the $`d_{z^2}`$ component of Mn $`e`$-orbital and the $`p_z`$ component of nearest neighbor N p orbital are the dominant states close to $`E_F`$. On the other hand, Mn $`t_2`$ and nearest neighbor As $`p`$ states hybridize strongly in (Mn,Ga)As. This is also true for (Mn,Ga)N in the zinc-blende structure. Symmetry of the crystal structure and the splitting of d-states under the corresponding crystal fields determine the position of $`t_2`$ and $`e`$ states. For the rest of the atoms in the unit cell, states having $`s`$ and $`p`$ character are the dominant.

## IV Interatomic exchange interactions

To determine the interatomic exchange interactions, we followed a simple model. In the unit cell, two Mn atoms were placed in various positions. For each Mn-Mn separation, ferromagnetic (FM) and antiferromagnetic (AFM) alignments of Mn spins were considered. The total energy difference $`\mathrm{\Delta }E`$ ($`\mathrm{\Delta }E`$=$`E_{tot}^{AFM}`$-$`E_{tot}^{FM}`$) between these two alignments is a measure of interatomic exchange interaction. In Fig. 6, we plot $`\mathrm{\Delta }E`$ as a function of Mn-Mn separation $`d`$ for both Mn doped GaAs and Mn doped GaN systems. For Mn doped GaN, both zinc-blende and wurtzite structures were considered for these calculations with two Mn atoms in unit cells having 64 atoms and 72 atoms respectively. For Mn doped GaAs, a 64 atom unit cell was chosen. Ferromagnetic interaction between Mn spins is favored for all the cases considered here. In (Mn,Ga)N, the first nearest neighbor (nn) exchange interaction is the strongest. The value of $`\mathrm{\Delta }E`$ is increased a little bit compared to the same for (Mn,Ga)As. So, for a defect-free calculation, there is no indication that (Mn,Ga)N should have a much higher $`T_C`$ than (Mn,Ga)As. It suggests that the formation of other phases during the growth is responsible for very high $`T_C`$ observed in certain experiments. It is also seen from the figure that $`\mathrm{\Delta }E`$ decreases sharply with $`d`$ for Mn doped GaN systems. It shows that ferromagnetic exchange interaction in (Mn,Ga)N is short-ranged. This indicates that the formation of Mn clusters within a short radial distance might lead to high values of $`T_C`$. For Mn doped GaAs, the exchange interaction is long ranged and doesn’t decrease rapidly. The proper range of ferromagnetic interactions can be studied with a bigger supercell. In summary, the results for exchange interactions indicate that the ferromagnetic interaction between Mn spins in (Mn,Ga)As is mediated by delocalized valence band holes whereas the origin of ferromagnetism in (Mn,Ga)N may result from a double-exchange mechanism sato  involving the hopping of Mn $`d`$ electrons. The other probable mechanism can be the formation of Zhang-Rice magnetic polaron zr .

Another interesting observation is the anisotropy of the exchange interactions in different crystallographic directions. In (Mn,Ga)As, ferromagnetic coupling is stronger either in the bonding direction e.g.  or in a direction where the two Mn spins are connected by As bonds e.g. in  direction. The coupling decreases for the  direction where there is no As atom in between (shown as a filled triangle in fig. 6). The strong ferromagnetic interactions in  or  directions result from strong $`pd`$ hybridization between Mn $`d`$ and As $`p`$ states. As the interaction is mediated by the delocalized states, they are sufficiently long-ranged. In  direction, the exchange interaction between the Mn spins can not be mediated by the $`pd`$ hybridization. So, the value of $`\mathrm{\Delta }E`$ is decreased compared to the other cases. In (Mn,Ga)N, the exchange interaction is not mediated by the delocalized valence band states and the ferromagnetic interactions decrease sharply with the Mn-Mn separation. So, the anisotropy in the exchange interactions is not significantly observed. We can conclude that in general, this anisotropy should be present for all Mn doped semiconductors where ferromagnetic long range interaction is mediated by delocalized valence band states. We also studied the exchange interactions in (Mn,Ga)N within the LSDA+U scheme in a similar way described above. U=4 eV and U=7 eV were considered for the calculations. We found that $`\mathrm{\Delta }E`$ for nearest neighbor Mn-Mn distance increases from 170 meV/Mn to 240 meV/Mn while going from U=4 eV to U=7 eV. Ferromagnetic interactions between Mn spins become stronger as the localization of the Mn $`d`$-states is increased.

## V Conclusion

We have studied the electronic structure and magnetism of Mn doped GaN systems for a wide concentration range of Mn. The deep acceptor level of Mn lies distinctly in the gap of GaN, separated from the valence band of host GaN. This is in contrast to the case of (Mn,Ga)As where the Mn forms shallow acceptor level close to the valence band of GaAs. Ferromagnetic interactions are short-ranged in (Mn,Ga)N systems whereas they have delocalized itinerant character for (Mn,Ga)As. Also, in (Mn,Ga)As, exchange interactions are significantly anisotropic in different crystallographic directions. This has been explained in terms of anisotropic $`pd`$ hybridization. The results presented here are for ideal systems having no contribution from the native defects formed during the non-equilibrium growth. To have a more realistic picture, one should take into account these effects. Also, Mn can occupy interstitial positions in the lattice and the spin interactions can become quite complicated. We are presently calculating the formation energies of the defects and subsequently the magnetic interactions in presence of them. Results will be reported in future communications.

###### Acknowledgements.

We are grateful to the Swedish Foundation for Strategic Research (SSF) , the Swedish Natural Science Research Council (NFR), and the Göran Gustafsson Foundation for financial support. We thank J. Hafner and G. Kresse for giving us the possibility to use VASP. B.S. is grateful to Dr. Lars Nordström and Prof. Olle Eriksson for valuable suggestions. Also, we acknowledge computational support from National Supercomputer Center (NSC) at Linköping, Sweden.