text:
# Synchronized Formation of Sub-Galactic Systems at Cosmological Reionization: Origin of Halo Globular Clusters

## 1. Introduction

Zinn (1985) pointed out insightfully that there appears to exist two separate populations of globular clusters in the Galaxy: one which is spherically distributed, has low metallicity with no radial metallicity gradient and low rotation, the other which is more concentrated towards the Galactic center and has high metallicity, a gradient in metallicity and high rotation. Such a dichotomy seems to be shared by other galaxies (e.g., Ashman & Bird 1993; Forbes, Grillmair, & Smith 1997; Forbes, Brodie, & Huchra 1997; Kundu 1999; Gebhardt & Kissler-Patig 1999; Barmby et al. 2000; Beasley et al. 2000). The first population, the halo globular clusters in both the Galaxy and external galaxies, display an additional array of interesting properties (Harris 1991). First, they are old with age $`10`$Gyr. Second, their characteristics are uniform across all galaxy types and sizes. Third, their mass function resembles a power-law ($`M^{1.7\mathrm{to}2.0}`$) at the high mass end.

It has been a long standing challenge to explain such features. Several important theoretical models including the pre-galactic Jeans mass based model (Peebles & Dicke 1968) and the secondary thermal instability based model (Fall & Rees 1985) have been put forth but may not be without some intriguing difficulties. For example, the Jeans mass based model is not naturally expected to produce power-law mass function for globular clusters, while the thermal instability based model appears to require delicate heating sources (Harris & Pudritz 1994; Meylan & Heggie 1997). A fully consistent model for the origin of the halo globular clusters remains unavailable at this time. In this paper, we suggest that the reionization of the universe could trigger a simultaneous formation of sub-galactic systems at that epoch, whose characteristics match remarkably well those of the observed halo globular clusters.

After describing the globular cluster formation theory in §2, we discuss some of the expected properties of globular clusters which form in the current theory in §3, followed by conclusions in §4. The following cosmological model is used throughout: $`\mathrm{\Omega }_{CDM}=0.26`$, $`\mathrm{\Omega }_b=0.04`$, $`\mathrm{\Omega }_t=\mathrm{\Omega }_{CDM}+\mathrm{\Omega }_b=0.3`$, $`\mathrm{\Lambda }=0.0`$, $`H_0=65`$km/sec/Mpc and $`\sigma _8=1.0`$; the global gas to total mass ratio $`R\mathrm{\Omega }_b/\mathrm{\Omega }_t=0.13`$.

## 2. Formation of Globular Clusters at Cosmological Reionization

In the standard picture of structure formation, the universe is thought to be reionized mostly by photons from quasars or stellar systems more massive than $`10^9\mathrm{M}_{}`$ (Haiman, Rees & Loeb 1997; Gnedin & Ostriker 1997). Less massive systems, after having produced a trace amount of stars (Pop III) at an earlier epoch, can no longer form stars due to lack of cooling processes (Haiman, Thoul, & Loeb 1996; Haiman, Rees & Loeb 1997; Tegmark et al. 1997). During the cosmological reionization phase small halos with neutral gas sitting idly within (i.e., without internal radiating sources) suddenly see a sea of radiation approaching. We will show that this radiation field has a dramatic dynamic effect on them. While an exact treatment of this process is quite intricate and would require detailed three-dimensional radiation-hydrodynamical simulations, we present a simple analysis that we believe captures the essence of the process. To make the calculations relatively tractable, we assume, without loss of primary characteristic features, that the halos are initially spherical with the following density profile:

$`\rho (r)`$ $`=\rho _v({\displaystyle \frac{r}{r_v}})^2\text{for}r<r_v`$ (1)
$`=\rho _v({\displaystyle \frac{r}{r_v}})^3\text{for}r>r_v,`$
where $`\rho _v100\rho _c(z)`$ is the density at the virial radius $`r_v`$ and $`\rho _c(z)=\rho _c(0)\mathrm{\Omega }_t(1+z)^3`$ is the critical density at redshift $`z`$ \[$`\rho _c(0)`$ is the critical density today\]. The virial mass of the halo (i.e., the total mass interior to $`r_v`$) is

$`M_v=4\pi \rho _vr_v^3`$ (2)
and the one-dimensional velocity dispersion within the virialized region ($`r<r_v`$) is

$`\sigma _v=\sqrt{2\pi G\rho _vr_v^2},`$ (3)
where $`G`$ is the gravitational constant. We assume that the ratio of gas to total mass in the halos under consideration is equal to the global baryon to total mass ratio $`R`$. The total mass within the radial range $`r_v`$ to $`r>r_v`$ (outside the virial radius) is

$`M(r_vr)=M_v\mathrm{ln}{\displaystyle \frac{r}{r_v}}.`$ (4)
The clumpy universe is reionized outside in: low density regions are ionized first and higher density regions (except those in the most massive, ionizing sources) become ionized at progressively later times (Miralda-Escudé, Haehnelt, & Rees 2000; Gnedin 2000b). The reionization phase progresses on a time scale of a Hubble time at the redshift in question; the mean (volume weighted) radiation field builds up slowly up to a value of approximately $`10^{24}`$erg/cm2/hz/sec/sr at Lyman limit. An accelerated phase follows, when about $`90\%`$ of the baryons have been ionized and a sudden jump in the amplitude of the mean radiation field intensity at Lyman limit to $`10^{22}10^{21}`$erg/cm2/hz/sec/sr occurs within a redshift interval of a fraction of unity and the reionization process is said to be complete (Miralda-Escudé et al. 2000; Gnedin 2000b). The mean free path of ionizing photons approaches the horizon size at a much later time. As Figure 3 (below) will show that the gas mass contained in the halos of interest here (including gas surrounding the halos that will be compressed) is of order $`50\%`$ of total gas, so we are concerned with the phase when a large fraction of the gas is not yet ionized and the radiation field is still increasing very slowly (see Figure 2a of Gnedin 2000b).

Let us now examine the global ionization front (I-front), which propagates from low density regions towards halos. The speed of the external global I-front towards a halo is

$`v_I={\displaystyle \frac{m_p}{f_HR}}J_\nu \rho _v^1\left({\displaystyle \frac{r}{r_v}}\right)^3`$ (5)
at $`r>r_v`$, where $`J_\nu _{\nu _0}^{\mathrm{}}4\pi j_\nu 𝑑\nu /h\nu `$ is the number of ionizing photons per cm2 sec and $`j_\nu `$ is the intensity of the meta-galactic radiation field; where $`f_H=0.76`$ is the hydrogen mass fraction of baryons; $`m_p`$ is the proton mass. Adopting a powerlaw $`j_\nu =j_{LL,21}10^{21}(\nu /\nu _0)^\beta `$erg/cm2/hz/sec/sr one obtains $`J_\nu =1.9\times 10^6\beta ^1j_{LL,21}`$cm-2 sec-1. Using $`\beta =2`$ and $`R=0.13`$ we obtain

$`v_I=1.8\left({\displaystyle \frac{r}{r_v}}\right)^3\left({\displaystyle \frac{1+z_{ri}}{16}}\right)^3\left({\displaystyle \frac{j_{LL,21}}{2\times 10^3}}\right)\text{km/s},`$ (6)
where $`z_{ri}`$ is the redshift when this occurs. Note that for a softer stellar spectrum of $`j_\nu `$, $`\beta 4`$ and $`v_I`$ would be lower.

The I-front is accompanied by a shock front. In order for the I-front to drive a strong compressive shock inward, the shock velocity has to exceed the velocity of the ionization front (D-type), which decreases with decreasing halo-centric radius. Otherwise, the neutral gas gets ionized and raised to a high pressure before the gas is reached by the shock front (see Spitzer 1978 for an introduction on the subject and Bertoldi & McKee 1990 for a detailed treatment in the context of interstellar clouds). We assume that this condition sets the initial radius of the inward shock, $`r_i`$: at $`r_i`$ the ionization front and the shock have the same velocity. In the present case the shock velocity will instantaneously settle to a speed of $`v_{shk}=4/3\sqrt{kT_{ri}/m_p}=12(T_{ri}/10^4K)^{1/2}`$km/s at the start of the shock due to the ram pressure (see equation 7 below), where $`T_{ri}`$ is the temperature of photo-ionized gas and $`k`$ is the Boltzmann constant. We will return to this when we discuss gas cooling in the shell later in this section. Since $`v_{shk}15`$km/s, this initial radius is outside the virial radius (see equation 6) (note that, for a more realistic density profile where the density slope near the viral radius is about $`2.4`$ instead of $`3`$ assumed here, the initial shock radius would be still larger).

We make some simplifying assumptions to treat such a shock. We assume that the shock is strong, which is justified since neutral gas has a very low temperature, and the material that is swept up along its way is accumulated in a thin shell obeying the momentum equation

$`{\displaystyle \frac{d(M_sv_s)}{dt}}={\displaystyle \frac{GM(<r)M_s}{r^2}}+4\pi r^2(p_{ext}p_{int}),`$ (7)
where $`M_s`$ is the gas mass within the shell; $`v_s`$ is the velocity of the shell; $`M(<r)`$ is the total mass interior to radius $`r`$; the first term and the two terms within the parentheses on the right hand side are due to gravity, external and internal pressure, respectively. In addition, we assume that dark matter does not respond to the motion of the gas, which is unlikely to be in serious error for our purpose. In any case, our assumption is conservative in the sense that inclusion of dark matter response to the gas collapse would enhance the collapse of halo gas.

It is instructive to qualitatively consider the calculation in two separate regions: $`r_vr<r_i`$ and $`0r<r_v`$. At $`r_v<r<r_i`$ we assume that the gas is initially cold (i.e., $`p_{int}`$ may be ignored) and $`p_{ext}=R\rho _v(r_i/r_v)^3kT_{ri}/m_p`$ is constant. The shell will be subject to external pressure and gravity, counter-balanced by ram pressure, and will obtain a certain velocity when reaching $`r=r_v`$. At $`0<r<r_v`$, if the shell travels supersonically, the downstream pressure is undisturbed. Since initially the gas in the virialized region at $`0<r<r_v`$ is in hydrostatic equilibrium, the region interior to the shell at any moment will remain in hydrostatic equilibrium. In other words, for a shell traveling supersonically, the gravitational force and interior pressure force just cancel out. In the absence of the external pressure and the self-gravity of the shell (which is small until the shell becomes self-gravitating), the shell will cruise inward (as long as it remains supersonic), except that it slows down due to the swept-up gas.

Quantitatively, we integrate the momentum equation (7) numerically from $`r_i`$ to $`r=0`$, with the initial condition $`v_s(r_i)=0`$ and $`M_s=0`$. We define $`\eta `$ in the following equation and require $`\eta `$ to be greater than $`1`$ (momentum condition) in order for the shell to reach the center of the halo:

$$\eta \frac{\mathrm{min}[v_s(r)]}{\sigma _v}1,$$
(8)
where $`\mathrm{min}[v_s(r)]`$ is the minimum velocity of the shell in the radial range $`[0,r_v]`$, and $`\sigma _v`$, the velocity dispersion, is used as the isothermal sound speed of the gas at $`r<r_v`$. Note that if equation (8) is not met and the shell slows down to a subsonic speed at some radius (likely at a large radius since most of the mass is at large radii), the interior gas will adjust its pressure causing the interior gas to be adiabatically compressed, which will quickly build up the internal pressure and cause the shell to slow down and eventually to stop.

Figure 1 shows $`\eta `$ (equation 8) as a function of $`\sigma _v`$ for $`T_{ri}=2\times 10^4`$Kelvin and $`j_{LL,21}=2\times 10^3`$ at $`z_{ri}=15`$. The adopted value of photoionized gas temperature is conservative \[Miralda-Escudé & Rees (1994) note that the post ionization temperature may be as high as $`5\times 10^4`$Kelvin\]. We see that halos with velocity dispersion $`\sigma _v11`$km/sec could collapse momentum-wise. The dashed curve is the total collapsed baryonic mass $`M_{gas}`$ (see the right vertical axis) and the dot-dashed curve is $`M_{GC}=0.2M_{gas}`$, where $`0.2`$ is a putative value for star formation efficiency (see below for further discussion). We see that compact, baryonic systems with mass $`M_{gas}10^6\mathrm{M}_{}`$ may form.

Let us ask how much compression is required to produce globular cluster like systems. Figure 2 shows the initial virial radius (solid curve) and the radius of the compressed gas cloud whose velocity dispersion is $`\sigma _{cm}=\sigma _v/0.2`$ (short dashed curve) as a function of halo velocity dispersion at $`z_{ri}=15`$. We see that the gas clouds need to be compressed by a factor of $`100`$.

A compressed gas cloud will be flattened (due to angular momentum) and shocked, cool, presumably fragment and form stars. However, we know little about high redshift star formation in such systems observationally or theoretically, and simply assume that the star formation efficiency is $`20\%`$. From an analysis of Local Group dwarf galaxies Gnedin (2000a) concludes that star formation at high redshift is about 4%, if it was continuous, but could be significantly higher, if most of the star formation occurred well before $`10`$Gyr ago. Therefore, our choice of star formation efficiency seems plausible but uncertainty may be as large as a factor of two to five. Note that the mean baryonic number density in the compressed cloud is of order $`1000`$cm-3, with the actual final density in a flattened disk probably even higher, matching that of molecular cloud cores (Friberg & Hjalmarson 1990; Goldsmith 1987). With such a large compression the dark matter would become a small fraction (of order $`1\%`$ for a more realistic halo density profile) of the final total mass.

The stars formed will likely blow away most of the remaining gas by supernova explosions (Dekel & Silk 1986; Mac Low & Ferrara 1999). After the remaining gas is lost as a wind, the cloud will expand by approximately a factor of $`1/0.2`$ and result in a final stellar system with a velocity dispersion of $`\sigma _{GC}=0.2\times \sigma _{cm}=\sigma _v`$ and a size shown as the long dashed curve in Figure 2. We see that the typical size of such a stellar system is between $`240`$pc. Both the assumption that the final stellar system has the same velocity dispersion as the initial host halo and the fudge factor $`1/0.2`$ in setting $`\sigma _{cm}`$ are rather uncertain, although these two parameters are degenerate and can be condensed into a single parameter. But small variations would not change the results qualitatively. Nevertheless, detailed simulations will be necessary to definitively quantify them.

It is yet unclear whether the required gas compression indicated in Figure 2 is achievable in real halos, even if the momentum condition $`\eta >1`$ (equation 8) is met. The critical issue here is that a finite angular momentum of the halo gas before compression (Peebles 1969; White 1984) could set a maximum compression factor. Let us denote the initial spin parameter of the halo (both gas and dark matter) as $`\lambda _iJ_i|E_i|^{1/2}G^1M_i^{5/2}`$ (where $`J_i`$, $`E_i`$ and $`M_i`$ are the initial angular momentum, total energy and total mass, respectively), then the final spin parameter of the gas cloud is (assuming no loss of angular momentum) $`\lambda _f\lambda _iR^1\sigma _{cm}/\sigma _v40\lambda _i`$ (for the short dashed curve in Figure 2). Therefore, a typical halo with $`\lambda _i0.05`$ (Barnes & Efstathiou 1987; Ueda et al. 1994; Steinmetz & Bartelmann 1995; Cole & Lacey 1996) would give $`\lambda _f2`$. Clearly, only halos with spin parameter values significant lower than the typical value can be adequately compressed to form the globular cluster like systems. Are there enough number of such low spin halos at $`z15`$?

Figure 3 shows the cumulative fraction of total mass in halos at $`z=15`$ (solid curve) using the Press-Schechter (1976) formula. The region left of the vertical (long dashed) line has $`\eta 1`$ (equation 8). We see that the (total) mass fraction of halos that meet the condition indicated by equation (8) is approximately $`15\%`$. As numerical simulations have shown that the distribution of the $`\lambda _i`$ is quite broad with the 80% range being $`0.010.10`$ (Ueda et al. 1994; Steinmetz & Bartelmann 1995; Cole & Lacey 1996); there are 10% of halos having $`\lambda _i0.01`$, which would give $`\lambda _f0.4`$ in this model \[when being compressed by a factor of $`100`$ to $`r_{cm}`$ (short dashed curve) in Figure 2\]. If $`10\%`$ of the halos form globular clusters, the resultant density parameter is $`\mathrm{\Omega }_{GC,comp}0.2\times 1.7\times \mathrm{\Omega }_b\times 15\%\times 10\%2\times 10^4`$ (where the five terms from left to right are the star formation efficiency, the ratio of the total swept-up gas mass over the gas mass within the virial radius, the baryon density, the mass fraction of halos in question and the fraction of low spin halos that collapse). The observed ratio of globular cluster mass over total Galactic baryonic mass today is approximately $`0.25\%`$ (Harris & Racine 1979), giving $`\mathrm{\Omega }_{GC,obs}=0.0025\mathrm{\Omega }_{}=1.5\times 10^5`$, where $`\mathrm{\Omega }_{}=0.006`$ is the stellar mass density (e.g., Gnedin & Ostriker 1992). Evidently, there are enough low spin halos in the indicated mass range that could account for the observed halo globular clusters. Since the slope of the predicted mass function is also in agreement with observations, one can infer that the expected number density of globular clusters is consistent with observations. The issue of exactly what the threshold of the spin parameter is in order to form a globular cluster system can only be settled by detailed simulations. Complications such as shell instability may play a role in this regard. It may be that the initial $`\mathrm{\Omega }_{GC}`$ is significantly larger than the presently observed value, because a large fraction (perhaps relatively faster spinning) of the initial globular clusters have largely evaporated or disrupted in time and now are a part of the halo stellar population.

The assumption of a thin shell needs further examination. With the usual shock jump conditions and equation (7), we find that at the start of the shock propagation the shell velocity settles to $`v_s=\sqrt{kT_{ri}/m_p}`$ and the shock front travels ahead of the shell with a velocity of of $`v_{shk}=4v_s/3`$. The postshock gas temperature would be $`T_s=1/3T_{ri}=6667`$Kelvin and the postshock pressure would be $`4p_{ext}/3`$. Therefore, while the shell is bounded at the inner side by the ram pressure, the outer boundary will slightly expand initially. In order to keep the shell thin the shell gas needs to be cooled. We compute the gas cooling including all relevant processes (Haiman et al. 1996) with the assumption that the gas has been subject to the radiation for a period of $`t_{Hubb}`$ (the Hubble time at $`z_{ri}`$) and the density is $`4\rho _v`$. For the gas in question the primary cooling process is molecular hydrogen cooling and the primary heating process is photo-heating of hydrogen atoms. Figure 4 shows the ratio of the cooling time to the time scale of the system in question $`t_{sys}`$ as a function of optical depth at the Lyman limit, where $`t_{sys}1\text{kpc}/10\text{km/sec}1\times 10^8`$yr. The incident radiation field has $`j_{LL,21}=2\times 10^3`$ with an index $`2`$ and the postshock gas temperature is assumed to be $`5\times 10^3`$Kelvin. We see that gas with a moderate optical depth (shielding) of $`\tau 410`$ can cool efficiently (gas at low end of the optical depth is heated up by photo-heating, while gas at the high end of the optical depth can not cool efficiently due to a low abundance of molecular hydrogen). For a shell of neutral gas we can relate its mass to the optical depth as $`M(\tau )=3.4\times 10^4\tau (r/1\text{kpc})^2\mathrm{M}_{}`$. Apparently, adequate shielding may be achieved for the gas and gas may be able to cool efficiently. Note that the gas cooling time remains approximately constant, assuming it cools isobarically, until it reaches a few hundred Kelvin. In addition, self-gravity of the shell will provide some confinement to keep the shell thin. More definitive answers can not be given without a detailed radiation-hydrodynamic simulation.

Finally, we give a dimensional analysis for the mass threshold of compressed gas clouds. We use the following line of reasoning. The upper mass threshold of the cloud should be related to the temperature of the external ionized gas by $`kT_{ri}=GM_vm_p/r_v`$, where the virial radius $`r_v`$ is related to $`M_v`$ by $`\delta \frac{4\pi }{3}r_v^3(1+z_{ri})\mathrm{\Omega }_t\rho _c(0)=M_v`$ ($`\delta 200`$ is the virial overdensity). The temperature of the external ionized gas should be related to hydrogen ionization potential $`kT_{ri}=\frac{1}{2}m_ec^2\alpha ^2`$, where $`m_e`$ is electron mass, $`c`$ is the speed of light and $`\alpha `$ is the fine structure constant. Combining these relations and the Friedmman equation we obtain $`M_{gas}=\frac{R}{2^{5/2}\pi }\frac{c^3\alpha ^3}{HG}(\frac{m_e}{m_p})^{3/2}=3\times 10^8R(\frac{1+z_{ri}}{15})^{3/2}\mathrm{\Omega }_t^{1/2}h^1\mathrm{M}_{}`$, where $`H`$ is the Hubble constant at $`z_{ri}`$ and $`R`$ is the global gas to total mass ratio. It is seen that the upper mass threshold is around $`10^{7.5}\mathrm{M}_{}`$, consistent with the detailed calculations given above.

## 3. Characteristics of the Globular Clusters

The dotted curve in Figure 3 is the computed mass function of halos, $`n(M)`$ (in units of $`h^3Mpc^3`$), using the Press-Schechter formula and the dashed curve a power-law mass function, $`n(M)dMM^2dM`$. A globular cluster mass function of $`n(M)dM=M^2dM`$ seems to be borne out naturally in this picture, in excellent agreement with observations for globular clusters more massive than $`10^5\mathrm{M}_{}`$ (Surdin 1979; Racine 1980; Richtler 1992). The observed shallower slope at lower mass may be due to disruptive effects during subsequent dynamical evolution of globular clusters, which tend to work against low mass systems (Chernoff, Kochanek, & Shapiro 1986; Chernoff & Shapiro 1987; Spitzer 1987; Aguilar, Hut, & Ostriker 1988; Chernoff & Weinberg 1990; Harris 1991; Goodman 1993; Gnedin & Ostriker 1997). We note, however, there may be a cutoff in the velocity dispersion of globular clusters at $`\sigma _{GC}2`$km/s, corresponding to $`\sigma _{cm}10`$km/s, below which cooling and thus star formation are less efficient at the absence of atomic cooling.

It is likely that a significant fraction of the initial globular clusters may possess a significant amount of angular momentum and their initial shapes may be significantly flattened. While much more work is needed to quantify the dynamical evolution of rotating many-body stellar systems (e.g., Einsel & Spurzem 1999), Akiyama & Sugimoto (1989) show that initially rotating stellar systems tend to become rounder with time. Although it is not possible at present to make useful comparisons with observations, some initial rotation of the globular clusters might be preferred by extant observations (e.g., White & Shawl 1987),

Simulations of Ueda et al. (1994) (also Cole & Lacey 1996) show a correlation between average spin parameter $`\overline{\lambda }_i`$ and halo mass, $`\overline{\lambda }_iM_v^{0.2\pm 0.1}`$, for massive (galaxy size) systems that they simulate. If we assume that for the halos in question a similar correlation exists, such as $`\overline{\lambda }_iM_v^{1/6}`$, and that gas clouds of all halos are compressed to a state of comparable final spin parameter, then the average final cloud radius $`\overline{r}_{GC}r_v\overline{\lambda }_i^2r_vM_v^{1/3}M_v^0`$ (because $`r_vM_v^{1/3}`$), i.e., the average final radius of the globular clusters is independent of the initial mass of the halo, in agreement with the observed lack of correlation between the half-light radius and the luminosity (Djorgovski & Meylan 1994). Of course, since there is a range in the initial spin parameter values, there will be a corresponding spread in the sizes of globular clusters, as observed. Since the final velocity dispersion of a globular cluster is $`\sigma _{GC}(M_v/\overline{r}_h)^{1/2}`$, the independence of $`\overline{r}_h`$ of $`M_v`$ dictates $`\sigma _{GC}M_v^{1/2}L^{1/2}`$, in agreement with the observed correlation $`\sigma _{GC}L^{0.6\pm 0.15}`$ (as well as the observed correlation with surface brightness $`\sigma _{GC}I_h^{0.45\pm 0.05}`$) (Djorgovski & Meylan 1994). A prediction that may be made is that the logarithmic dispersion in $`\sigma _{GC}`$ in the $`\sigma _{GC}L`$ correlation should be about half the size of the logarithmic dispersion in the size of globular clusters, which seems to be in agreement with observations (Djorgovski & Meylan 1994).

There appears to exist a tendency for halos with relatively large $`\lambda `$ to reside preferentially in low density regions in simulations (e.g., Ueda et al. 1994). This effect would conceivably produce the observed correlation of the specific frequency of globular clusters with galaxy type in that early type galaxies have a high number of globular clusters per unit galactic luminosity than late type galaxies (Harris 1991). This arises because early type galaxies originate from initial higher density peaks and would have (on average) lower spins and thus higher specific frequencies. Quite intriguing is that this trend might have some bearing on the formation of spiral and elliptical galaxies and Hubble sequence. But within a similar galaxy type, the number of globular clusters in a galaxy may simply be proportional to the mass or luminosity of the galaxy, as observed (Hanes 1977; Harris & Racine 1979), since the mass of a galaxy basically indicates the size of the initial (comoving) region from which it has collected matter hence proportionally the number of globular clusters. More certain is that the range of properties of globular clusters such as size, luminosity and velocity dispersion should be quite uniform across all galaxies, independent of size, type, age, etc, as observed (Harris 1991). This is due to the fact that the formation of globular clusters in the present model has little correlation with later large-scale density fluctuations that form galaxies (except the correlation between spin and overdensity mentioned above).

Chemical enrichment by Pop III stars, which is computed to be $`10^{3.5}\mathrm{Z}_{}`$ (Ostriker & Gnedin 1996), is much lower than the observed stellar metallicity of metal poor halo globular clusters of $`10^{1.2}\mathrm{Z}_{}`$ (Zinn 1985). It thus appears that some low level of self-enrichment has to occur in order to account for the observed metallicity of globular clusters. While a detailed treatment of star formation in the proposed globular clusters is not possible, the two-generation star formation scenario (Cayrel 1986; Brown et al. 1991,1995; Zhang & Ma 1993; Parmentier et al. 1999) could work in this picture. We note that the escape velocity of a compressed cloud (where star formation occurs) is larger than its corresponding globular cluster by a factor of $`5`$ (this is rather uncertain, though). Therefore, confinement of the first generation of supernovae (Dopita & Smith 1986) may be more readily accommodated than in other scenarios. Then, the second generation of stars formed in a massive starburst (with the assumed formation efficiency of 20%) completely blows away the remaining gas, except for perhaps the most massive globular clusters (Dekel & Silk 1986; Mac Low & Ferrara 1999), for which further self-enrichment may take place. This is in accord with the observed internal chemical homogeneity of most globular clusters (cf. Larson 1988), with the exception of the most massive systems such as $`\omega `$ Cen where some chemical inhomogeneities are observed (Cohen 1981).

While halo globular clusters depicted here should be, on average, metal poor, as observed (Zinn 1985), complex star formation histories quite likely introduce a large dispersion in metallicity among globular clusters. Thus it may be expected that larger galaxies that formed by collecting matter in a larger (comoving) region are expected to display a larger variance in metallicity, in agreement with observations (Harris 1991). Since the metallicity distribution is asymmetric and bounded by $`Z=0`$ at the low end, a large variance gives a larger mean, interestingly consistent with observations (Harris 1991). On the other hand, since metallicity does not play a noticeable role in the formation of individual globular clusters (i.e., dynamic collapse of the gas clouds) in this model, it should be expected to show little correlation with any other quantities of globular clusters, such as luminosity, size, velocity dispersion, in agreement with observations (Djorgovski & Meylan 1994).

If the universe went through a second, separate reionization phase of singly ionized helium at lower redshift ($`z3`$), as recent observations hint (e.g., Reimers et al. 1997), a similar radiation I-front will sweep through halos of neutral gas. However, due to much lower density and lower abundance of helium (compared to hydrogen) at lower redshift, the speed of the I-front is likely to be larger than the shock speed until well within the virial radius; i.e, the I-front will remain as an R-type and no significant compression will occur. Therefore, our model indicates that the age spread of metal poor globular clusters should be less than $`1`$Gyr, consistent with recent observations (Rosenberg et al. 1999).

Finally, it seems inevitable that a fraction of the original globular clusters may be in the intergalactic space where no galaxies were formed to collect them. These intergalactic globular clusters are possibly at a somewhat younger dynamic state than Galactic halo globular clusters due to lack of external dynamical effects. The exact number density of intergalactic globular clusters depends sensitively on the angular momentum distribution of halos in low density regions and requires more detailed calculations for quantification. No doubt it will be important to search for them with the next generation of large telescopes.

## 4. Conclusions

Radiation from stars in galaxies with mass $`M_t>10^9\mathrm{M}_{}`$ or quasars are thought to be mostly likely responsible for reionizing the universe at redshift $`z620`$. However, a large fraction of mass resides in density concentrations on smaller mass scales, i.e., gas rich sub-galactic halos with mass $`M_t<10^8\mathrm{M}_{}`$. These smaller systems can not cool to form stars without an external trigger. Reionization of the universe has an interesting dynamic effect on these halos. The rise of the external radiation field causes a synchronous propagation of an inward, convergent shock towards each halo. We show that, for halos with a velocity dispersion $`\sigma _v11`$km/s, the resident gas of mass $`M_b10^410^7\mathrm{M}_{}`$ will be compressed. For a significant fraction of these halos that have relatively low initial angular momentum ($`\lambda 0.01`$), a compression of a factor of $`100`$ in radius is possible and compact self-gravitating baryonic systems would form. Such gaseous systems fragment to form stars and become stellar systems of size $`240`$pc, velocity dispersion $`110`$km/s and a total stellar mass of $`M_s10^310^6\mathrm{M}_{}`$.

Their expected properties all seem to be in agreement with the observed salient features of halo globular clusters, as discussed in the previous section. Angular momentum distribution of small halos at high redshift plays a critical role in determining many of the properties of globular clusters. The expected number density of halo globular clusters is consistent with what is observed. The mass function of slope $`2`$ at the high mass end is predicted by the model. Strong correlation between velocity dispersion and luminosity (or surface brightness) and lack of correlation between velocity dispersion and size are expected. Globular clusters should display similar properties, regardless of the host galaxy type, age, size or luminosity. Metallicity is, on average, expected to be low and should not correlate with any other quantities of globular clusters, except that a larger dispersion of metallicity among globular clusters is expected for larger galaxies. A trend of specific frequency with galaxy type may be produced in the model. We propose that these stellar systems represent the initial state of the presently observed halo globular clusters.

It will be important to systematically search for high redshift as well as intergalactic globular clusters, because they will provide tests of the proposed model and could potentially shed light on the theory of universe reionization and the general formation picture of galaxies. Also important is to make detailed radiation-hydrodynamic simulations with high resolution and superior shock resolving power to study the reionization process and follow the dynamics of gas collapse. Both are formidable challenges.

I thank Bruce Draine, Jerry Ostriker and Michael Strauss for a critical reading of the manuscript and many cogent comments, Stu Shapiro for suggesting a dimensional analysis, Zoltan Haiman for allowing me to use his cooling program, and Rennan Barkana and Jordi Miralda-Escudé for useful discussion. This research is supported in part by grants AST93-18185 and ASC97-40300.


format: text

text:
# Tunneling magnetoresistance in diluted magnetic semiconductor tunnel junctions

## I INTRODUCTION

Since the early experiment done by Julliere, there has been considerable interest in studying the tunneling magnetoresistance (TMR) effect in magnetic tunnel junctions  due to its promising applications in digital storage and magnetic sensor technology. The conventional electrodes of magnetic tunnel junctions are made up of the ferromagnetic transition metals such as Fe, Co, Ni, and their alloys. Usually the TMR effect in the conventional magnetic tunnel junctions is described by the Julliere model. The electrons near at the Fermi levels of the electrodes participate in the tunneling process and have the main contribution to the tunneling conductance and so to the TMR effect.

Recent discovery of the ferromagnetism in the diluted magnetic semiconductors (DMS) has opened the rich field of an interesting interplay between magnetic and semiconducting properties, which makes it possible to combine the information processing and data storage in one material. At high concentration $`c`$ of randomly distributed $`\mathrm{Mn}^{2+}`$ ions doped in GaAs samples with high hole density $`c^{}`$, Ga1-xMnxAs compounds exhibit ferromagnetism with transition temperature as high as 110 K at certain value of $`x`$. The origin of ferromagnetism in the DMS can be explained by the itinerant holes through the $`p`$-$`d`$ exchange coupling between the valance-band holes and the spins of impurity Mn ions. In order to investigate the spin-polarized electron transport phenomenon in the DMS ferromagnetic semiconductor heterostructures, the epitaxial tunnel junctions have been successfully fabricated in the form of Ga0.961Mn0.039As/AlAs/Ga0.961Mn0.039As. The TMR effect was observed and the maximum TMR ratio thus obtained is as high as 44.4% at 4.2 K. The subsequent experiment  also realized the TMR effect in the DMS tunnel junctions of Ga0.95Mn0.05As/AlyGa1-yAs/Ga0.97Mn0.03As.

The Julliere model  with several necessary extensions  is known to describe the tunnel junctions, in which the two electrodes are made of the conventional ferromagnetic metals or the half-metallic manganites such as La1-xSrxMnO3. For the manganite tunnel junctions, various tunneling processes are included in the model as follows. The elastic spin-flip scattering of the spin polarized carriers by the local magnetic impurities within the insulating barrier and the inelastic spin-flip scattering induced by the thermal collective excitations of local spins at the interfaces between the insulating barrier and the manganite electrodes are the two basic ingredients. By taking into account the temperature dependence of the spin polarization of itinerant carriers, one can calculate the temperature dependence of the TMR ratio. In the paper, we have applied the above generalized Julliere model to the DMS tunnel junctions and obtained the various properties of the TMR effects of the system.

The paper is organized as follows. In Sec. II, we propose a spin-polarized tunneling model for the DMS tunnel junctions by generalizing the previous spin-polarized tunneling model in order to correctly take into account the basic physics of the DMS. In Sec. III, we show that the experimentally observed TMR ratio in the DMS tunnel junctions can be well reproduced from our calculations with the appropriate choice of typical material parameters and also show that the TMR ratio thus obtained is strongly dependent on both the itinerant-carrier density and the magnetic ion density in the DMS electrodes. Finally we conclude with a brief summary in Sec.IV.

## II SPIN-POLARIZED TUNNELING MODEL

We will briefly summarize the basic physics of ferromagnetism occurring in certain DMS compounds following the seminal work of König et al..

In the doping range of our interest, that is, $`0.035<x<0.06`$, Ga1-xMnxAs can be described as a degenerate semiconductor with relatively high hole concentrations and the $`\mathrm{Mn}^{2+}`$ ions act as local magnetic impurities.

The Hamiltonian describing the basic physics of the DMS can be written as

$$=d^3r\underset{\sigma }{}\mathrm{\Psi }_\sigma ^{}(𝐫)(\frac{\mathrm{}^2}{2m^{}}^2\mu ^{})\mathrm{\Psi }_\sigma (𝐫)$$

$$+J_{pd}d^3r𝐒(𝐫)𝐬(𝐫),$$
(1)
where $`𝐒(𝐫)`$ stands for the spin density of the local magnetic impurities, which is given by $`_I𝐒_I\delta (𝐫𝐑_I)`$ with $`𝐑_I`$ denoting the positions of magnetic impurities. The field operator $`\mathrm{\Psi }_\sigma ^{}(𝐫)`$ ($`\mathrm{\Psi }_\sigma (𝐫)`$) is the creation (annihilation) operator for the itinerant hole at position $`𝐫`$ with spin $`\sigma `$. The $`𝐬(𝐫)`$ represents the spin density of itinerant-carriers, which can be written in terms of field operators as $`\frac{1}{2}_{\sigma \sigma ^{}}\mathrm{\Psi }_\sigma ^{}(𝐫)\stackrel{}{\tau }_{\sigma \sigma ^{}}\mathrm{\Psi }_\sigma ^{}(𝐫)`$, where $`\stackrel{}{\tau }`$ are the three Pauli spin matrices. The $`m^{}`$ and $`\mu ^{}`$ are the effective mass and chemical potential of itinerant carriers, respectively. The $`J_{pd}`$ stands for the $`p`$-$`d`$ exchange coupling strength between itinerant holes and $`\mathrm{Mn}^{2+}`$ ion impurity spins. The recent experiments suggest that the itinerant hole density $`c^{}`$ is much lower than the magnetic ion density $`c`$ satisfying $`c^{}/c1`$. Here we focus our attention on the DMS with ferromagnetic ground state. Using the Holstein-Primakoff transformations and performing the coarse graining procedures, the local spin density $`𝐒(𝐫)`$ can be written in terms of the newly defined bosonic fields $`b^{}(𝐫),b(𝐫)`$ as follows

$$S^+(𝐫)=\sqrt{2cSb^{}(𝐫)b(𝐫)}b(𝐫),$$
(2)
$$S^{}(𝐫)=b^{}(𝐫)\sqrt{2cSb^{}(𝐫)b(𝐫)},$$
(3)
$$S^z(𝐫)=cSb^{}(𝐫)b(𝐫).$$
(4)
The self-consistent spin wave approximation developed in Ref.  yields the following effective Hamiltonian for the DMS system

$$_{\mathrm{eff}}=\underset{𝐤\sigma }{}\epsilon _{𝐤\sigma }c_{𝐤\sigma }^{}c_{𝐤\sigma }+\underset{|𝐪|<q_c}{}\mathrm{\Omega }_𝐪b_𝐪^{}b_𝐪,$$
(5)
where $`\epsilon _{𝐤\sigma }`$ stands for the band energy of the itinerant holes with spin $`\sigma =\pm 1`$, which is given by

$$\epsilon _{𝐤\sigma }=\epsilon _𝐤+\frac{\mathrm{\Delta }}{2}\sigma $$
(6)
and $`\epsilon _𝐤=\mathrm{}^2k^2/(2m^{})`$. The spin wave dispersion $`\mathrm{\Omega }_𝐪`$ is given by

$$\mathrm{\Omega }_𝐪=\frac{J_{pd}}{2}(n_{}n_{})\frac{J_{pd}\mathrm{\Delta }}{2V}\underset{𝐤}{}[f(\epsilon _𝐤)f(\epsilon _{𝐤+𝐪})]$$

$$\times \frac{1}{\mathrm{\Omega }_𝐪+\epsilon _{𝐤+𝐪}\epsilon _𝐤},$$
(7)
where $`c_𝐤`$ and $`b_𝐪`$ are the itinerant hole and the spin wave annihilation operators in the momentum space respectively, and $`f(x)`$ is the Fermi-Dirac distribution function. Here $`n_{}`$ ($`n_{}`$) is the spin-down (spin-up) itinerant hole density satisfying $`n_{}+n_{}=c^{}`$, the exchange gap $`\mathrm{\Delta }=cJ_{pd}S`$, and $`q_c`$ is the Debye cutoff for the spin waves of local impurity spins with $`q_c^3=6\pi ^2c`$. At finite temperatures, the spin wave dispersion can be generalized by imposing the following self-consistency conditions for the finite temperature exchange gap : $`\mathrm{\Delta }(T)=J_{pd}S^z`$. Here $`S^z`$ represents the thermal average of the impurity Mn ion spins, which are approximately calculated by the following formula

$$S^z=\frac{1}{V}\underset{|𝐪|<q_c}{}S_s(\beta S\mathrm{\Omega }_𝐪)$$

$$=\frac{1}{V}\underset{|𝐪|<q_c}{}\{Sn_B(\mathrm{\Omega }_𝐪)+(2S+1)n_B[(2S+1)\mathrm{\Omega }_𝐪]\},$$
(8)
where $`_s(x)`$ is the Brillouin function, $`n_B(x)`$ the Bose-Einstein distribution function, and $`\beta =1/k_BT`$. The second term in the second equality of Eq. (8) describes how the thermally induced spin waves from the independent Bose statistics with no limit in the number of spin waves indeed reduce the magnetization of the system. The third term takes into account the correct spin kinematics, which rules out the unphysical states.

Based on the above scheme, one can calculate the magnetization of the Mn impurity spins and the band structure of the itinerant holes approximately self-consistently by solving a set of coupled equations. In this theoretical framework, the non-monotonic dependence of the Curie temperature on the itinerant hole density $`c^{}`$ has been successfully obtained.

Now we consider the DMS magnetic tunnel junctions in which both electrodes are made of the ferromagnetic DMS Ga1-xMnxAs compounds. The carrier tunneling through a thin AlAs insulating barrier involves removing a carrier at one side and creating it at the other side with spin conserving or spin-flipping. The spin-flip tunneling processes can be divided into two different categories: The elastic spin-flip tunneling arising from the $`\mathrm{Mn}^{2+}`$ ion impurity spins located within the insulating barrier and the inelastic spin-flip tunneling induced by the thermal collective excitations of the impurity spins near at the interfaces between the insulating barrier and the DMS electrodes.

For the collinear magnetization configurations, the tunneling Hamiltonian of the DMS tunnel junctions can be written as

$$H=H_\mathrm{L}+H_\mathrm{R}+H_\mathrm{T},$$
(9)
$$H_\mathrm{L}=\underset{𝐩\sigma }{}\epsilon _{𝐩\sigma }c_{𝐩\sigma }^{}c_{𝐩\sigma }+\underset{|𝐪|<q_c}{}\mathrm{\Omega }_𝐪b_{𝐪L}^{}b_{𝐪,L},$$
(10)
$$H_\mathrm{R}=\underset{𝐤\sigma }{}\epsilon _{𝐤\sigma }d_{𝐤\sigma }^{}d_{𝐤\sigma }+\underset{|𝐪|<q_c}{}\mathrm{\Omega }_𝐪b_{𝐪R}^{}b_{𝐪R},$$
(11)
$$H_\mathrm{T}=\underset{\mathrm{𝐤𝐩}\sigma }{}(T_{\mathrm{𝐤𝐩}}^0d_{𝐤\sigma }^{}c_{𝐩\sigma }+iT_{\mathrm{𝐤𝐩}}^{}d_{𝐤\sigma }^{}c_{𝐩,\sigma }+\mathrm{H}.\mathrm{c}.)$$

$$+\frac{1}{\sqrt{V}}\underset{\mathrm{𝐤𝐩},|𝐪|<q_c}{}\{T_{\mathrm{𝐤𝐩𝐪}}^{\prime \prime }[S_L^z(𝐪)+S_R^z(𝐪)](d_𝐤^{}c_𝐩d_𝐤^{}c_𝐩)$$

$$+\mathrm{H}.\mathrm{c}.\}+\frac{1}{\sqrt{V}}\underset{\mathrm{𝐤𝐩},|𝐪|<q_c}{}\{T_{\mathrm{𝐤𝐩𝐪}}^{\prime \prime }[\sqrt{2cS_L}b_{𝐪L}+\sqrt{2cS_R}b_{𝐪R}]$$

$$\times (d_𝐤^{}c_𝐩+c_𝐩^{}d_𝐤)+\mathrm{H}.\mathrm{c}.\}.$$
(12)
Here $`H_\mathrm{L}`$ ($`H_\mathrm{R}`$) is the effective Hamiltonian for the bulk DMS in the left (right) electrode, which has been derived from the original Hamiltonian of Eq. (1) as explained earlier. The operator $`d_{𝐤\sigma }`$ $`(c_{𝐩\sigma })`$ is the carrier annihilation operator with momentum $`𝐤(𝐩)`$ and spin $`\sigma `$ of the right (left) DMS electrode, and $`S_\alpha ^z(𝐪)=cS_\alpha b_{𝐪\alpha }^{}b_{𝐪\alpha }`$ with $`\alpha =L`$, $`R`$.

The term with coefficient $`T_{\mathrm{𝐤𝐩}}^{}`$ in $`H_\mathrm{T}`$ represents the spin-flip tunneling arising from the Mn2+ ion impurity spins inside the AlAs insulating barrier. This type of spin-flip tunneling has been introduced to the tunneling Hamiltonian for the conventional ferromagnetic tunnel junctions and the manganite tunnel junctions.  The imaginary $`i`$ for the spin-flip tunneling matrix elements indicates that the elastic spin-flip tunneling processes are incoherent in any magnetization configurations. As a generalization of the typical tunneling Hamiltonian  to the DMS tunnel junctions, the last two terms have been included in the tunneling Hamiltonian of Eq. (12), which describe the inelastic spin-flip tunneling process induced by the thermal spin wave collective excitations of local spins near at the interfaces between the insulating barrier and the DMS electrodes. Since this kind of carrier tunneling involves the emission or absorption of the local spin collective excitations near the left or right interfaces, it is an inelastic tunneling process. The inelastic spin-flip tunneling may occur not only directly at the interfaces but also within the coherence length of the interfaces, which is relatively long comparing to the barrier width. It then allows us to use the bulk spin wave dispersion relation in the tunneling Hamiltonian. The fact that the wave vector $`𝐤`$ is independent of $`𝐩`$ in Eq. (12) implies that we only consider the incoherent tunneling processes, which are known to be appropriate for the defect-populated insulating barriers and the rough interfaces. As usual, we make the following assumptions that the strengths of various tunneling amplitudes, that is, $`|T_{\mathrm{𝐤𝐩}}^0|^2,|T_{\mathrm{𝐤𝐩}}^{}|^2`$, and $`|T_{\mathrm{𝐤𝐩𝐪}}^{\prime \prime }|^2`$ are independent of momenta $`𝐤,𝐩,𝐪`$ and then will be taken to be their average values $`|T^0|^2,|T^{}|^2`$, and $`|T^{\prime \prime }|^2`$, respectively. By making a transformation for the tunneling Hamiltonian from the collinear magnetization configurations to arbitrary relative angle $`\theta `$ between the spin orientations of the two electrodes and using the standard Green’s function methods, one can obtain the tunneling conductance $`G(\theta )`$ at zero bias up to the leading order in the tunneling amplitudes, which is given by

$$G(\theta )=\frac{\pi e^2}{\mathrm{}}|T^0|^2[\rho _{}(\epsilon _F)+\rho _{}(\epsilon _F)]^2[1+P^2\mathrm{cos}\theta $$

$$+\gamma (1P^2\mathrm{cos}\theta )],$$
(13)
where $`P`$ stands for the itinerant carrier spin polarization given by

$$P=\frac{\rho _{}(\epsilon _F)\rho _{}(\epsilon _F)}{\rho _{}(\epsilon _F)+\rho _{}(\epsilon _F)},$$
(14)
and $`\rho _{}(\epsilon _F)`$ ($`\rho _{}(\epsilon _F)`$) is the density of states (DOS) of the spin-down (spin-up) carriers at the Fermi level of the DMS electrodes. The relative angle $`\theta `$ may be tunable by applying a magnetic field and $`\gamma `$ is given by the following relation

$$\gamma =\gamma _1+4\eta cS(cSS^z),$$
(15)
where $`\gamma _1=|T^{}|^2/|\stackrel{~}{T}|^2`$ and $`\eta =|T^{\prime \prime }|^2/|\stackrel{~}{T}|^2`$ with $`|\stackrel{~}{T}|^2=|T^0|^2+2c^2S^2|T^{\prime \prime }|^2|T^0|^2`$, since $`c^2|T^{\prime \prime }|^2`$ is usually about two orders of magnitude smaller than $`|T^0|^2.`$ In deriving Eq. (15), we have used the following relation $`V^1_{|𝐪|<q_c}b_𝐪^{}b_𝐪=cSS^z`$. The corresponding tunneling resistance is given by $`R(\theta )=1/G(\theta )`$.

The TMR ratio is defined by the following formula

$$\left(\frac{\mathrm{\Delta }R}{R_\mathrm{P}}\right)_{\mathrm{max}}=\frac{R_{\mathrm{AP}}R_\mathrm{P}}{R_\mathrm{P}},$$
(16)
where $`R_{\mathrm{AP}}`$ and $`R_\mathrm{P}`$ represent the tunneling resistance in the antiparallel and parallel spin alignments of the two ferromagnetic electrodes respectively, and the maximum TMR ratio for the DMS tunnel junctions can be obtained as

$$\left(\frac{\mathrm{\Delta }R}{R_\mathrm{P}}\right)_{\mathrm{max}}=\frac{2(1\gamma )P^2}{1P^2+\gamma (1+P^2)}.$$
(17)
If the electrodes of the DMS tunnel junctions are made of the non-degenerate DMSs, the tunneling conductance will be proportional to the itinerant carrier density and the TMR ratio will be more appropriately given by Eq. (17) with the carrier spin polarization replaced by $`P=(n_{}n_{})/(n_{}+n_{})`$. The different definitions of the carrier spin polarization in the TMR ratio for the degenerate and non-degenerate DMS tunnel junctions are due to the distinctly different transport properties in the corresponding DMSs.

## III results and discussions

In Fig. 1 and Fig. 2 are illustrated the temperature dependences of the impurity Mn ion magnetization of the DMS, the itinerant-carrier spin polarization, and the maximum TMR ratio for the DMS tunnel junctions at several different itinerant-carrier densities $`c^{}`$ with the following choice of material parameters: $`\gamma _1=0.692`$, $`\eta =0.005\mathrm{nm}^6`$, and $`m^{}=0.5m_e`$, $`J_{pd}=0.15\mathrm{eV}\mathrm{nm}^3`$ appropriate for the bulk Ga1-xMnxAs compounds  and fixed impurity Mn ion density $`c=1\mathrm{nm}^3`$, where $`m_e`$ is the free electron mass. Figure 1 (a) shows the impurity Mn ion magnetization $`S^z`$ normalized by $`c`$. One can notice that the critical temperature $`T_c`$ decreases with the further increase of the itinerant carrier density $`c^{}`$ from $`c^{}=0.1\mathrm{nm}^3`$. From the solid curves of Fig. 1 (b) corresponding to $`c^{}=0.1\mathrm{nm}^3`$, one can see the following features that at relatively low temperatures, the spin-up carrier band is completely lifted from the spin-down carrier band by the large exchange gap $`\mathrm{\Delta }(T)`$. At the mean-field level, this leads to the complete spin polarization of the itinerant carriers, i.e. $`P=100\%`$, exhibiting the half-metallic feature. In this case the TMR ratio is reduced to $`(\mathrm{\Delta }R/R_p)_{\mathrm{max}}=1/\gamma 1`$, which decreases with the increase in the population of the thermally excited spin wave collective excitations leading to the increase of $`\gamma `$ as shown in Fig. 1 (c). As the temperature increases, the size of the exchange gap $`\mathrm{\Delta }(T)`$ is gradually reduced. With the further increase of temperature beyond a certain value, the number of the majority spin-down carriers begins to decrease, while that of the minority spin-up carriers starts to increase. At high temperatures, the carrier spin polarization decreases rapidly with temperatures, which makes the TMR effect disappear above a critical temperature. The factor of $`(1\gamma )`$ in Eq. (17) indicates that the TMR effect may vanish below the critical temperature, which is due to the competition between the spin-conserving tunneling amplitude $`|T^0|^2`$ and the spin-flip tunneling ones $`|T^{}|^2,|T^{\prime \prime }|^2`$.

For higher carrier densities, say, $`c^{}=0.2\mathrm{nm}^3`$, the Fermi energy of the system increases and then both spin states are occupied losing the half-metallic feature as shown in Fig. 1 (b). Since the partial carrier spin polarization is not beneficial to the TMR effect, it leads to the low TMR ratio. Figure 2 (a) illustrates the impurity Mn ion magnetization $`S^z`$ normalized by $`c`$ at different values of $`c^{}`$. In contrary to Fig. 1 (a), one can see that $`T_c`$ decreases with the decrease of the itinerant carrier density $`c^{}`$ from $`c^{}=0.05\mathrm{nm}^3`$ exhibiting the non-monotonic dependence of $`T_c`$ on $`c^{}`$, which has been emphasized in Ref. . For $`c=1\mathrm{nm}^3`$, the optimal carrier densities, which will achieve the maximum TMR ratio, are shown to be in the range of $`0.05\mathrm{nm}^3<c^{}<0.1\mathrm{nm}^3`$.

Recently Hayashi et al. have performed an experiment for Ga0.961Mn0.039As/AlAs/Ga0.961Mn0.039As tunnel junctions and obtained the TMR ratio of about 44% at 4.2 K. By fitting the above experimental data to our calculations, we were able to estimate the approximate values of material parameters: $`c^{}=0.13\mathrm{nm}^3`$, $`c=0.86\mathrm{nm}^3`$, and $`\gamma _1=0.692`$. The above value of carrier density $`c^{}`$ is about 15 $`\%`$ of the Mn ion concentration $`c`$. The value of $`\gamma _1`$ indicates that there exists strong spin-flip scattering by the Mn ions within the insulating barrier. During the fabrication process of the DMS tunnel junctions, the highly doped Mn ions in the Ga1-xMnxAs layers can easily diffuse into the AlAs insulating layer. Reducing the Mn ions in the insulating layer is crucial to achieve a large TMR ratio in the DMS tunnel junctions.

Figure 3 shows the temperature dependences of the impurity Mn ion magnetization of the DMS, the itinerant-carrier spin polarization, and the maximum TMR ratio for the DMS tunnel junctions at several values of impurity Mn ion densities $`c`$ with fixed carrier density $`c^{}=0.1\mathrm{nm}^3`$. Figures 3 (a) and (b) demonstrate that the carrier band structure is mainly determined from the size of the exchange gap $`\mathrm{\Delta }(T)`$. At low impurity densities, say, $`c=0.6\mathrm{nm}^3`$, the strength of the ferromagnetic coupling is weak, which induces partial carrier spin polarization. By increasing $`c`$, the ferromagnetic coupling becomes strong, which completely polarizes the carrier spins above a certain value of $`c`$. So achieving the high doping of the impurity Mn ions is necessary to have much higher $`T_c`$ and larger TMR ratio at high temperatures.

When the magnetic ion density is beyond a certain large value, say, $`c=1.5\mathrm{nm}^3`$, the impurity ion magnetization and the itinerant-carrier spin polarization seem to exhibit the behaviors of the first order ferromagnetic transition as shown in Fig. 3 (a) and (b). In this case, the first order ferromagnetic transition is induced by the $`p`$-$`d`$ exchange coupling, in contrast to that for the two-dimensional quantum well system at very low carrier densities, which is argued to be due to the carrier-carrier correlation effects.

Our results clearly demonstrate that the TMR effects in the DMS tunnel junctions are dominated by both the ferromagnetic coupling strength and the carrier band structure through the impurity ion density and the itinerant carrier density in the DMS electrodes. The elastic and inelastic spin-flip tunnelings play an important role in reducing the TMR ratio.

## IV summary

In summary, we have presented the TMR formula for the DMS tunnel junctions by generalizing the previous spin-polarized tunneling model and taking into account the basic physics of the DMS. The experimentally observed TMR ratio was reproduced by the appropriate choice of the typical material parameters. It was also shown that the TMR ratio has a strong dependence on the itinerant-carrier density and the magnetic ion density in the DMS electrodes as well. In view of the recent predictions of the room temperature ferromagnetic semiconductors, the DMS tunnel junctions can be practically very useful in the future digital storage and magnetic sensor technology. Our analysis of the TMR ratio can be potentially useful to achieve larger TMR ratio by optimally adjusting the material parameters.

## V acknowledgments

This work was supported by the Brain Korea 21 Project and also by Grant No. 1999-2-112-001-5 from the interdisciplinary Research program of the KOSEF.

Fig. 1. (a) Impurity Mn ion magnetization of the DMS, (b) itinerant-carrier spin polarization $`P`$, and (c) the maximum TMR ratio in the DMS tunnel junctions as a function of temperature for several different itinerant-carrier densities $`c^{}=0.3,0.2,0.18,0.15,0.1\mathrm{nm}^3`$ with the following choice of material parameters $`m^{}=0.5m_e`$, $`J_{pd}=0.15\mathrm{eV}\mathrm{nm}^3`$, $`\gamma _1=0.692`$ and $`\eta =0.005\mathrm{nm}^6`$, and the fixed value of $`c=1\mathrm{nm}^3`$.

Fig. 2. (a) Impurity Mn ion magnetization of the DMS, (b) itinerant-carrier spin polarization $`P`$, and (c) the maximum TMR ratio in the DMS tunnel junctions as a function of temperature at several different itinerant-carrier densities $`c^{}=0.05,0.03,0.01,0.001\mathrm{nm}^3`$ with the following choice of material parameters $`m^{}=0.5m_e`$, $`J_{pd}=0.15\mathrm{eV}\mathrm{nm}^3`$, $`\gamma _1=0.692`$ and $`\eta =0.005\mathrm{nm}^6`$, and the fixed value of $`c=1\mathrm{nm}^3`$.

Fig. 3. (a) Impurity Mn ion magnetization of the DMS, (b) itinerant-carrier spin polarization $`P`$, and (c) the maximum TMR ratio in the DMS tunnel junctions as a function of temperature for several different impurity Mn ion densities $`c=0.6`$, 0.7, 0.9, 1.2, 1.5 $`\mathrm{nm}^3`$ with the following choice of material parameters $`m^{}=0.5m_e`$, $`J_{pd}=0.15\mathrm{eV}\mathrm{nm}^3`$, $`\gamma _1=0.692`$ and $`\eta =0.005\mathrm{nm}^6`$, and the fixed value of $`c^{}=0.1\mathrm{nm}^3`$.