format: text

text:
# The evolution of primordial black hole masses in the radiation-dominated era

## I Introduction

The search for PBHs is certainly one of the most interesting areas of high-energy astrophysics. While quite strong limits have been published \[LG\] on the fractional mass in PBHs, it is possible that a fraction of short GRBs may correspond to evaporating holes having masses at formation of the order the Hawking mass $`M10^{15}g`$ that would be evaporating today \[CG\]. The celebrated Hawking radiation would be driving black holes to their end, although several important details regarding this final phase are still being analyzed. However, Sivaram \[Siv\] has argued that these objects would not be observed at all, since that these PBHs would have absorbed a substantial amount of mass when the Universe was very young and hot, and therefore would not be evaporating today. This cautionary remark is important, since then short GRBs must be produced by other sources. Here we reanalyze the problem, and show that substantial mass gain did not occur at all (Section 2). The kinematic analysis of Section 3 indicates that a possible boost to enhance the mass absorption term is severely suppressed. The last Section is devoted to a demonstration of the impossibility of keeping thermodynamic equilibrium between the black holes and the radiation, thus justifying the previous discussion and the expected evolution of the mass. A slightly more general argument along these lines is given in the Appendix.

## II  Evaporation and absorption of mass

Hawking \[SH\], found in a seminal paper that black holes would loss their mass through a slow but irreversible quantum process, where virtual particles from the vacuum are converted in real ones created by the large gravitational energy. This quantum process gives rise to mass loss where the black hole had its temperature increased continuously. The mathematical analysis of this problem, showed that this emitted spectrum is thermal, with an associated temperature,

$$T_{bh}=\frac{\mathrm{}c^3}{8\pi Gk_BM}\frac{10^7K}{(M/M_{})}.$$
(1)
Using the Stefan-Boltzmann law and the formula above, the evolution of the mass of the PBH is given by

$$\frac{dM}{dt}=\frac{A(M)}{M^2};$$
(2)
where $`A(M)`$ counts the degrees of freedom of the emitted particles. Numerically, $`A(M)`$ is given by

$$A(M)10^{24}g^3s^1\eta (M);$$
(3)
with $`\eta (M)=7.8\times D_1(M)+3\times D_{1/2}(M)`$, $`D_{1/2}`$ and $`D_1`$ are the multiplicity of quantum numbers (spin, charge and color) of the emitted particles. For the standard model with three fermionic generations, we have $`D_1=27`$ and $`D_{1/2}=90`$ at $`T100GeV`$ and $`A(M)7.8\times 10^{26}g^3s^1`$ (ignoring contributions from Higgs scalars and gravitons, see \[KT\] for details).

We know that the cosmological environment was very hot and dense in the radiation era, therefore, we expect some classical absorption of the energy-matter from the surroundings of a given PBH. We adopt a Friedmann-Robertson-Walker model, with $`R(t)`$ the scale factor describing the expansion (or contraction) of the universe and $`K=1,0,+1`$ are the possible values of the curvature parameter, describing, respectively open models ($`K=1`$ or $`K=0`$) and closed models ($`K=1`$). The dynamics of the FRW model is governed by the well known equations

$$\left(\frac{\dot{R}}{R}\right)^2=\frac{8\pi G}{3}\varrho \frac{K}{R^2};$$
(4)
$$\frac{\ddot{R}}{R}=\frac{4\pi G}{3}(\varrho +3P);$$
(5)
without cosmological constant (which is not relevant for our discussion). In the radiation-dominated era the universe was close to $`K=0`$ and the scale factor and the radiation temperature evolve as $`R(t)=R_i\times t^{1/2}`$ and $`T_{rad}(t)R(t)^1`$ respectively. Numerically, the temperature is given by $`T_{rad}1MeV(t/1s)^{1/2}`$ and then the radiation density is given by

$$\varrho _{rad}(t)8.4\times 10^4(t/1s)^2gcm^3.$$
(6)
We define as usual the cosmological horizon by the integral $`R_{hor}(t)=R(t)lim_{t_i0}_{t_i}^t\frac{cdt^{}}{R(t^{})}`$. The cosmological horizon exists when this limit converges for finite $`t`$. Taking the limit $`t_i0`$ we have for the solution $`R(t)t^{1/2}`$

$$R_{hor}(t)=2ct.$$
(7)
The mass contained in the horizon is given simply by the product $`M_{hor}(t)\frac{4\pi }{3}R_{hor}(t)^3\varrho (t)`$. Using eqs.(6) and (7) we obtain

$$M_{hor}(t)7.6\times 10^{37}(t/1s)g.$$
(8)
This relation is very important to the problem, since we expect that PBHs would be formed with masses comparable to the horizon mass \[Carr\]. There is also a dependence with the adiabatic index of the background matter $`\mathrm{\Gamma }=\frac{P}{\varrho }`$, an important feature for the collapse of perturbations associated to PBHs and generally overlooked. Since PBHs with initial masses larger than the horizon mass are excluded by causality, for a generic PBH with initial mass $`M`$ was formed in the early universe its formation time must satisfy

$$t_f(M)\mathrm{\hspace{0.17em}3}\times 10^{38}(M/g)s.$$
(9)
Therefore, from $`t_f(M)`$, we expect that the PBH mass will be below the horizon mass at subsequent times $`t>t_f`$.

In order to address the complete problem of the evolution of PBHs masses, we shall include a classical absorption term. This term must be included to the r.h.s. of eq.(2) and is of form

$$\left(\frac{dM}{dt}\right)_{abs}=\sigma _g(M)F_{rad}(T);$$
(10)
where $`\sigma _g(M)=\frac{27\pi }{4}r_{g}^{}{}_{}{}^{2}`$ is the gravitational cross section for the capture of relativistic particles \[ZN\]. Since $`F_{rad}(T)=c\varrho _{rad}(T)`$ is the radiative flux as seen in the rest frame of the black hole, the complete differential equation for the mass is given by

$$\frac{dM}{dt}=\frac{A(M)}{M^2}+\frac{27\pi G^2}{c^3}\varrho _{rad}(T)M^2;$$
(11)
where we have assumed that $`M>M_{Planck}`$ in order to avoid complicated corrections due to quantum gravity effects and neglected the back-reaction (see \[Norma\] and references therein for a quantum mechanical evaluation).

The inclusion of the second term is a necessary ingredient for the thermodynamical description, since PBHs formed obeying the causality constraints with masses $`M_{Haw}`$ before $`t1s`$ will experience an inward heat flux from the hotter surroundings. Consequently, these cold PBHs will gain mass-energy at the expense of the radiation. When the physical temperature of these PBHs equals to the radiation temperature, these objects will be in instantaneous equilibrium and $`\dot{M}=0`$ at this instant. We can derive the mass for equilibrium to happen at any epoch solving $`\dot{M}=0`$ \[NOS\]

$$M_c(t)\frac{10^{26}g}{(T/T_0)};$$
(12)
where $`T_0`$ is the present cosmic temperature. Any PBH with mass greater than this value at $`t`$ must be accreting ($`\dot{M}(M>M_c)>0`$) and not evaporating, since the Hawking evaporation is negligible in this case.

To be quantitative we may consider the ratio of the absorption to evaporation terms to obtain

$$\xi (M)=\frac{\dot{M}_{abs}}{\dot{M}_{evap}}=\left(\frac{M}{M_c}\right)^4;$$
(13)
which shows that the Hawking term quickly becomes negligible for masses above $`M_c`$.

A first look at eq.(13) could convince us to accept the conclusions of Sivaram \[Siv\], namely that these supercritical ($`M>M_c`$) PBHs would be growing very fast in the early universe. To show that this would be misleading we shall solve the eq.(11) for initially supercritical PBHs, (that is, ignoring the first term) in the thermal radiation bath evolving as stated in eq.(6). Then, the differential equation to solve is

$$\frac{dM}{dt}=\frac{27\pi G^2}{c^3}\varrho _{rad}(t)M^2.$$
(14)
The integration of this equation is immediate and we obtain

$$\frac{1}{M(t)}=\frac{1}{M_i}+1.2\times 10^{39}\left(\frac{1}{t}\frac{1}{t_i}\right)g^1;$$
(15)
where $`M_i`$ is the initial mass (at formation) and $`t_i`$ is the formation time (measured in seconds) for this object. This solution is valid if the initial mass satisfies $`M_c(t_i)<M_iM_{hor}(t_i)`$, as stated above.

We may invert the analytical solution eq.(15) above, casting it in the form

$$M(t>>t_i)=\frac{M_i}{[1ϵ(M_i,t_i)]};$$
(16)
with $`ϵ(M_i,t_i)=1.2\times 10^{39}\frac{(M_i/g)}{(t_i/1s)}`$.

Since $`ϵ(M_i,t_i)<<1`$ for PBHs satisfying $`M_i<<M_{hor}(t_i)`$; and even in the most favourable circumstances $`ϵ(M_i,t_i)M_{hor}(t_i)`$, we can write eq.(16) expanding it in a Taylor series and find that the growth of these objects is $`MM_i[1+ϵ(M_i,t_i)]`$. Note that the mass gain is larger for larger masses, and also increases for early formation times for a given mass. In fractional terms

$$\frac{\delta M}{M}=\frac{MM_i}{M_i}=ϵ(M_i,t_i).$$
(17)
Note that the fractional mass gain saturates with time for $`t>>t_i`$, as expected from the time dependence of the radiation density $`\varrho _{rad}(t)t^2`$. Even at $`t=t_i`$ the universe was not dense enough to drive a very strong growing. We can show it in the following way. Suppose that the hot environment can indeed support a relativistic PBH mass gain of the gravitational radius $`\dot{r_g}c`$. Let us show that the actual mass accretion was many orders of magnitude below this value, i.e. $`\dot{M}<<\frac{c^3}{2G}2\times 10^{38}gs^1`$. For this purpose, assume that initially, we had $`\dot{M}=\frac{27\pi G^2}{c^3}\varrho _{rad}M^2\frac{c^3}{2G}`$. Solving it, we find the minimal density of the radiation capable to drive the relativistic growing for that mass $`M`$

$$\varrho _{min}5.7\times 10^{51}gcm^3(M/M_{Haw})^2.$$
(18)
If we compare $`\varrho _{min}`$ with eq.(6) we obtain that these PBH would have to be formed before $`t_{rel}(M)3.8\times 10^{24}(M/M_{Haw})s`$. However, according to eq.(9) causal PBH formation is allowed after $`t_f>3\times 10^{23}(M/M_{Haw})s`$. Thus, most of the time the actual growing has to be much smaller than the maximal theoretical value $`\dot{M}_{max}=\frac{c^3}{2G}`$, which needs very extreme conditions.

A brief relativistic gain may nevertheless be allowed, since $`t_f`$ and $`t_{rel}`$ are very close. This means that a PBH with horizon mass, formed exactly at $`t_f`$ will gain some mass, but the maximal fractional gain is limited to be $`0.04`$. Therefore, if we consider the maximal gain, a PBH formed with $`M_i10^{15}g`$ would finish its evolution with final mass (at the end of the radiation era) greater than the Hawking mass, but with a small gain. However, that is enough to delay its evaporation, since the evaporation timescale is given by

$$t_{evap}(M)t_0(M/M_{Haw})^3;$$
(19)
therefore the mass $`MM_i+0.04M_i`$ then we obtain $`t_{evap}(M)(1.04)^3t_0`$.

The situation is quite different for PBHs with initial mass of order $`M_i`$ smaller than $`M_{Haw}`$. In this case, even the maximal gain can not drive these PBHs above the Hawking mass at all. Therefore, these objects will show very small modifications of their behaviour. Note that the maximum gain is easily obtained for $`M_iM_{hor}`$. Then substituting we have for that absolute maximal gain

$$ϵ_{max}(M_iM_{hor})1.17\times 10^{39}M_i\left[\frac{1}{3\times 10^{38}M_i}\right]0.04.$$
(20)
In order to understand the full behaviour of the PBHs masses it is useful to consider the evolution of the critical mass. The latter may be expressed in the form

$$M_c(t)10M_{Haw}(t/1s)^{1/2}.$$
(21)
If we compare the behaviour of a PBH with initial mass above the critical mass (but obviously below the horizon mass) and the growth of the critical mass, we observe that this object will cross the critical mass some time after its formation. While this PBH stays above the critical mass, it will not evaporate, but its final mass (that is, before the crossing time) will be given by eq.(19). The crossing time $`t_{cross}`$ is defined by

$$M(t)=M_i[1+ϵ(M_i,t_i)]10M_{Haw}(t_{cross}/1s)^{1/2};$$
(22)
with solution

$$t_{cross}(M_i)10^2s(M_i/M_{Haw})^2[1+ϵ(M_i,t_i)]^2.$$
(23)
This explains why these objects did not gain much mass: the time interval in which these objects stay above the critical mass is very small if compared to cosmological times, or to the timescale for evaporation $`t_{evap}t_0M^3`$.

Any supercritical PBH will be driven by the second term in eq.(14) and any subcritical PBH will be driven by the first term, the Hawking evaporation. In the section 3 we will show that the objects crossing the critical mass curve can not stay there, and must drift to the evaporation region of Fig.1 in which the evaporation quickly dominates.

## III Mass gain induced by high proper motions

Since we expect that at least for large Lorentz factors the absorption and evaporation processes must be modified, we want to consider corrections due to the proper motion of the PBH relative to the background radiation. This effect is described by the following expression

$$T_{rad}(\theta ,v)=\frac{T_{rad}(\theta ,v=0)}{\gamma (1\stackrel{}{n}.\stackrel{}{\beta })};$$
(24)
where we will denote $`T_{rad}(\theta ,v=0)`$ simply as $`T_{rad}`$, the temperature of the background radiation in the rest frame of the PBH when this PBH has $`v=0`$. Then, $`T_{rad}(\theta ,v)`$ for arbitrary $`v`$, denotes the angular distribution of the background radiation seen from the rest frame of this PBH when it moves with peculiar velocity given by $`v`$. The angle $`\theta `$ is defined in the scalar product $`\stackrel{}{n}.\stackrel{}{\beta }=v\beta cos(\theta )`$. $`\stackrel{}{n}`$ is the vector defining the sight of view for an arbitrary observer located at the Schwarzschild radius, $`\gamma =\frac{1}{\sqrt{1\beta ^2}}`$ and $`\beta =(v/c)`$.

The radiation density has an angular distribution as seen from a stationary observer located at the black hole horizon. The sign in the denominator of eq.(24) is such that we have a blueshift in the direction of the motion and a redshift in the opposite. The absorbed radiation is an angular average and it is given by \[NOS\]

$$<\varrho _{rad}(T_{rad})>=\frac{1}{4\pi }𝑑\mathrm{\Omega }\varrho _{rad}(T_{rad}(\theta ,\varphi )).$$
(25)
This integral is performed to give

$$\left(4\gamma ^2\frac{1}{\gamma }\right)\varrho _{rad}(T_{rad}(t));$$
(26)
where the time dependence in $`T_{rad}(t)`$ arises from the cosmological expansion, i.e. $`T_{rad}R^1`$.

Inserting this result in the eq.(11), we arrive at the differential equation for the mass of a black hole moving with velocity $`v`$, as seen by an asymptotic observer. Evidently $`\varrho _{rad}`$ in eq.(26) must be the same quantity given by eq.(6). Now, we seek a $`\gamma `$ capable to drive the rate of mass accretion to the maximal rate $`\dot{M}_\gamma \frac{c^3}{2G}`$. Scaling the black hole mass to Hawking mass units

$$\left(4\gamma ^2\frac{1}{\gamma }\right)=\frac{7\times 10^{53}}{(M/M_{Haw})^2}(\varrho _{rad}(t)/gcm^3)^1.$$
(27)
This is a second degree equation for $`\gamma `$ with solution

$$\gamma (t)=\frac{\alpha }{8\mu ^2\varrho _{rad}(t)}\left[1+\sqrt{1+\frac{16}{\alpha ^2}M^4\varrho _{rad}(t)^2}\right];$$
(28)
where $`\alpha =7\times 10^{53}`$, $`\mu M/M_{Haw}`$ and $`\varrho _{rad}(t)`$ is expressed in $`gcm^3`$. Substituting $`\varrho _{rad}(t)=\varrho (t_f)(t/t_f)^2`$, we can rewrite eq.(28) as

$$\gamma (t)=\frac{\alpha (t/t_f)^2}{8\mu ^2\varrho _{rad}(t_f)}\left[1+\sqrt{1+\frac{16}{\alpha ^2}\mu ^4\varrho _{rad}(t_f)^2(t/t_f)^4}\right];$$
(29)
since $`\varrho _{rad}(t)8.4\times 10^4(t/1s)^2`$, and our solution holds only for $`t>t_f(M)`$ we obtain the following inequality

$$\varrho _{rad}(tt_f)4.4\times \frac{10^{51}}{\mu ^2}.$$
(30)
Therefore, the prefactor in eq.(28) has to be greater than $`18`$, for any $`tt_f`$. We may rewrite the required Lorentz factor as

$$\gamma (t)>18\left[1+\sqrt{1+\frac{16}{\alpha ^2}\mu ^4\varrho _{rad}(t)^2}\right].$$
(31)
Applying the same arguments to second term, we can finally write

$$\gamma (t)=18(t/t_f)^2\left[1+\sqrt{1+6\times 10^4(t/t_f)^4}\right].$$
(32)
For asymptotic times $`t>t_f`$, we have $`\gamma (t)36(t/t_f)^2`$.

Therefore, we conclude that even at the beginning of its life ($`tt_f`$) and when the PBH was formed with $`MM_{hor}`$, a high value for the proper velocity $`v`$ would be needed to drive the PBH into the fast growing mode (differing from $`c`$ in one part per thousand). The additional terms in the square brackets are positive and enhance the value of the minimum $`\gamma `$ in eq.(31). The physical reason is very simple: from the initial time formation, the radiation density decreases with time, and consequently we will need to impart very high values for the proper velocity to drive an object above the critical mass. This particular situation was certainly unlikely in the early universe, and we do not expect to have PBHs with substantial accretion of mass. It is worth to remark that the ingredient that allowed us to perform these numerical evaluations and conclusions is the concept of critical mass (see also \[NOS2\] for a related discussion).

## IV Impossibility of thermodynamical equilibrium between PBHs and the background radiation

In this section we address the problem whether black holes and the surrounding radiation can be in thermodynamical equilibrium. This would imply that the PBH mass to be constant. Page \[Dany\] analyzed this problem in a gedanken experiment with one black hole contained in a finite box plus thermal radiation at a constant volume. It is easy to understand that this black hole within this box will settle down in stationary equilibrium, with its mass remaining constant in time. In this case, the total energy of the radiation $`+`$ black hole system is

$$E=AT^4V_{box}+M_{bh}.$$
(33)
Since that this system is closed, we have for variations in mass and temperature

$$\delta M=4AV_{box}T^3\delta T.$$
(34)
If the black hole loses mass due to Hawking effect $`\delta M<0`$, then the radiation temperature must raise $`\delta T_{rad}>0`$ and drive the black hole into equilibrium. If the black hole absorbs too much mass the radiation temperature diminishes substantially and additional accretion will not be allowed. We see immediatly that the situation is more complicated if the box has variable volume. It is immediate to conclude that microscopic variations in the critical mass must obey

$$\delta M_c(T_{rad})^5\frac{\delta M}{V_{box}}.$$
(35)
For very high temperatures, the variations in the critical mass are very small, as expected.

We shall now analyze the more complex and realistic situation, consisting in black holes and radiation within a cosmological model. For this, we emphasize that the critical mass is a cosmological quantity, defined by the radiation temperature. Numerically, we have

$$M_c(t)\frac{10^{26}g}{(T_{rad}(t)/T_0)}.$$
(36)
Therefore, since the temperature decreases with cosmological time, the critical mass grows monotonically. However, we note that the critical mass allows an instantaneous equilibrium, since that $`\dot{M}(M=M_c)=0`$. As the universe expands a black hole initially in equilibrium with the background radiation must deviate from it because of the different temporal behaviours.

To prove this statement, let us evaluate the variation of the area of the black hole. We know that $`R_g=\frac{2GM}{c^2}`$, then the horizon area is given by $`A=4\pi R_{g}^{}{}_{}{}^{2}`$, and therefore we have

$$dA=\frac{32\pi }{c^4}G^2MdM.$$
(37)
The sign of the mass variation is related to the energy that crossed the event horizon as $`dE=c^2dM`$. Then an energy flux through the event horizon induces a variation in the area given by

$$\frac{dE}{dA}=\frac{c^6}{32\pi G^2M}.$$
(38)
Then, the time derivative of the energy flux per unit area will be given by

$$\left(\frac{d^2E}{dtdA}\right)=\frac{c^6}{32\pi G^2M^2}\frac{dM}{dt}.$$
(39)
Substituting eq.(11) into eq.(39), we obtain

$$\left(\frac{d^2E}{dtdA}\right)=\frac{c^6}{32\pi G^2}\left[\frac{A}{M^4}\frac{27\pi G^2}{c^3}\varrho _{rad}(t)\right].$$
(40)
We will denote eq.(40) as $`\sigma (M,t)`$ for simplicity. We see that the energy is zero when $`M=M_c`$, then $`\sigma (M=M_c,t)=0`$ for any $`t`$. However, by evaluating the time derivative at $`M=M_c`$ and we will prove that the first derivative at the critical mass is not zero, indicating an inevitable deviation from the initial equilibrium.

Note that $`\sigma (M,t)`$ is positive-defined for hot black holes (subcritical black holes) and negative in the opposite case. Now, the time derivative of $`\sigma (M,t)`$ at the critical mass is given by

$$\dot{\sigma }(M=M_c,t)=\frac{27}{32}c^3\dot{\varrho }_{rad}(t)>0;$$
(41)
and it is positive-defined since that the radiation density decreases with time. If we integrate eq.(40) with the initial condition $`M=M_c(t_i)`$ until $`t>t_i`$ and given that the radiation diluted during this interval due to expansion, we obtain the result

$$_{\sigma (M,t)}^{\sigma (M_c,t_i)}𝑑\sigma =\frac{c^6}{32\pi G^2}_{t_i}^t𝑑t\left[\frac{4A}{M^5}\frac{dM}{dt}+\frac{27\pi G^2}{c^3}\dot{\varrho }_{rad}(t)\right];$$
(42)
and the right side can be rearranged in the form $`\frac{c^6A}{8\pi G^2}_{M_c}^M\frac{dM}{M^5}\frac{27}{32}c^3_{\varrho _{rad}(t_i)}^{\varrho _{rad}(t)}𝑑\varrho _{rad}`$. Then we arrive at

$$\sigma (M,t)=\sigma (M=M_c,t_i)+\frac{Ac^6}{32\pi G^2}\left[\frac{1}{M^4}\frac{1}{M_{c}^{}{}_{}{}^{4}}\right]\frac{27}{32}c^3\mathrm{\Delta }\varrho _{rad}.$$
(43)
The first term $`\sigma (M=M_c,t_i)=0`$ by definition, and the last term is negative as we already know. If we consider that $`M=M_c+\delta M`$, the term in brackets can be evaluated to give

$$\left[\frac{1}{M^4}\frac{1}{M_{c}^{}{}_{}{}^{4}}\right]=\frac{1}{M_{c}^{}{}_{}{}^{4}}\left[1+\frac{1}{1+F(M_c,\delta M_c)}\right];$$
(44)
where $`F(M_c,\delta M_c)=(\frac{\delta M}{M_c})^4+6(\frac{\delta M}{M_c})^2+4(\frac{\delta M}{M_c})^3+4(\frac{\delta M}{M_c})`$. Therefore, the term in brackets is always negative, for arbitrary variations in $`\delta M`$, after the PBH had crossed the critical mass at $`t_i`$. Note that $`F`$ is negative because the terms containing higher powers of $`\frac{\delta M}{M_c}`$ are positive and they dominate the other terms (which are negative, since $`\delta M`$ is arbitrary). For large variations of $`\delta M`$ of any sign, we have $`F<0`$, then we deduce that $`\sigma (M,t)`$ is positive, and thus conclude that this PBH migrates to the evaporating regime.

Now, even if the variation $`\delta M`$ is assumed to be small and negative, we could imagine that $`F(M_c,\delta M)`$ may change its signal and become positive. We shall prove that this is never the case. If the net variation in $`M`$ is small, we can ignore the first term in the r.h.s. of eq.(40) and integrate to give

$$\sigma (M,t_i+\mathrm{\Delta }t)+\left(\frac{c}{t_i}\right)^3\varrho _{rad}(t_i)\delta t;$$
(45)
where $`\delta t`$ is a small time interval satisfying $`\delta t\frac{27\pi G^2c^3}{g}(M)\mathrm{\Delta }\varrho _{rad}`$ with $`g(M)=\left[4A^2/M^7AB\varrho _{rad}\right]^1`$.

In order to show that a PBH in an initial equilibrium $`M=M_c`$ deviates from it is enough to show that the first derivative $`\dot{\sigma }`$ is non-zero at $`M=M_c`$ and is always positive for arbitrary variations $`\delta M`$. We proceed by steps.

Since that we are interested in the algebraic sign of $`\sigma (M,t)`$ only, we change variables to absorb the coefficient $`\frac{c^6}{32\pi G}`$ in eq.(40) and rewrite $`\sigma `$ as

$$\sigma (M,t)=\frac{A}{M^4}B\varrho _{rad}(t).$$
(46)
To prove that $`\dot{\sigma }(M=M_c,t)>0`$ it is enough to see that its derivative is just $`\dot{\sigma }(M,t)=\frac{A}{M^5}\dot{M}B\dot{\varrho }_{rad}(t)`$. However, the first term is null at $`M=M_c`$ and the second term is positive because we assumed that the universe expands, then $`\dot{\sigma }(M,t)>0`$ beyond the critical mass. Suppose that the variation in $`M`$ is such that we write $`M=M_c+\delta M`$ and the ”final” mass is below the critical mass by a significant amount, say $`M_0`$. Therefore, the energy flux conserves its original sign from the initial equilibrium. Now, let us prove that the sign of $`\dot{\sigma (M,t)}`$ does not change if the small variation $`\delta M`$ is negative. For this purpose, we write $`\dot{\sigma }(M,t)=\frac{A}{M^5}\left(\frac{A}{M^2}+B\varrho _{rad}M^2\right)B\dot{\varrho }_{rad}(t)`$. If we substitute $`M=M_c+\delta M`$ taking into account that this change is very small, we can expand this expression and arrive at

$$\dot{\sigma }(M,t)=\frac{A}{M^5}\left[\frac{A}{M_{c}^{}{}_{}{}^{3}}+2M_cB\varrho _{rad}(t)\right]\delta M>0.$$
(47)
Therefore, the rate of change of $`\sigma (M,t)`$ drives the hole to the evaporation. To see that this condition of eq.(47) is enough we may calculate $`\frac{d\sigma }{dt}(M_c+\delta M)>0`$, thus implying that $`d\sigma (M_c+\delta M)>0`$ or $`\sigma (M_c+\delta M)\sigma (M_c)>0`$. This in turn means that $`\sigma (M_c+\delta M)>0`$, since $`\sigma (M_c)=0`$ and $`\sigma (M)>0`$ happen in the evaporation region, or when the PBH losses much more energy than it is able to absorb from the environment. In the opposite case, if the mass variation is large but negative, we have $`\sigma (M=M_c+\delta M)=\frac{A}{M^4}B\varrho _{rad}(t)\frac{A}{M^4}>0`$, and the first term dominates because $`\xi <<1`$ in this case. This completes the analysis for all negative variations $`\delta M`$ around the critical mass.

Now, we have in mind the general behaviour of the PBH mass evolution as the universe expands. From the eq.(46) we deduce that

$$\dot{\sigma }(M,t)=\frac{4A}{M^5}\dot{M}B\dot{\varrho }_{rad}(t).$$
(48)
For subcritical black holes, $`\dot{M}<0`$ and we know that $`\dot{\varrho }_{rad}(t)<0`$ due to cosmological expansion. Then, for subcritical PBHs is immediate that $`\dot{\sigma }>0`$.

In the opposite case of supercritical black holes, the first term is $`\frac{4A}{M^5}\dot{M}<0`$. Then, the algebraic signal of $`\dot{\sigma }(M,t)`$ is determined by the difference $`B\dot{\varrho }_{rad}(t)\frac{4A}{M^5}\dot{M}`$. However, we found in section 1 that the relative mass gain is strictly limited to a few percent. This means that the second term is negligible and the algebraic signal of $`\dot{\sigma }(M,t)`$ remains positive.

The algebraic sign of $`\dot{\sigma }(M,t)`$ is always positive in ”natural” circumstances. Some paths where $`\dot{\sigma }(M,t)<0`$ are not a priori excluded, but these paths correspond to very special conditions where $`\dot{\varrho }_{rad}>0`$ or the negative variation in $`\varrho _{rad}(t)`$ is very small. For realistic models of the radiative era, we know that $`\varrho _{rad}(t)t^2`$, then, this variation is not small compared to PBH mass gain. All those paths where $`\dot{\sigma }(M,t)<0`$ are entropically disfavoured.

In the general case when the universe expands, the natural evolution for a PBH in the radiative era is to follow the type of (almost) horizontal path drawn in the Figure 1, before crossing the critical mass curve. There is a small mass gain along these paths, as we explained previously. The paths are thermodynamically favourable and terminate at the instantaneous equilibrium point $`M_c(t=t_c)`$ at the crossing time $`t_c(M_i)`$ (marked with a cross). Now, since that $`\dot{\sigma (M,t)}>0`$ at the critical mass, the deviation of equilibrium is unavoidable. The PBHs follow the type downward path afterwards, where $`\sigma (M,t)>0`$ and its derivative are positive as already discussed. In other words, the Hawking evaporation will eventually dominate the classical absorption.

In the general case we should analyze the full behaviour of

$$\sigma (M(t),tt_c)=\sigma (M_i,t_i)+\frac{Ac^6}{32\pi G^2}\left[\frac{1}{M^4(t)}\frac{1}{M_{i}^{}{}_{}{}^{4}}\right]\frac{27c^3}{32}\mathrm{\Delta }\varrho _{rad}(t);$$
(49)
where $`M(t)`$ is a solution of $`\dot{M}=\frac{A}{M^2}+B\varrho _{rad}(t)M^2`$. If the initial conditions are given by $`M_i>M_c(t_i)`$, $`\varrho _{rad}(t)t^2`$ (and it is implicit that $`\varrho _{pbh}(t)<<\varrho _{rad}(t)`$ as before); and substituting the appropriate solution for $`M(t)`$ into eq.(49) above, we deduce afte taking $`t>>t_c(M_i)`$ that $`\sigma (M(t>>t_c),t)>0`$. Therefore such a black hole crossed the critical mass curve and is now below the critical mass at $`t`$. The hole is then hotter than the environment and must evaporate.

All these thermodynamical considerations show quite generally that an explosive growth of PBHs in the radiation era is forbidden. This conclusion in unavoidable when the problem is analyzed using carefully the critical mass. The validity of Sivaram’s solution would also lead to the conclusion that an initially sub-horizon PBH would have to cross the horizon in a short time given by $`\mathrm{\Delta }t\frac{(M_{hor}M_{pbh})}{3\times 10^{38}g}s`$. This kind of behaviour was already ruled out by a set of calculations \[ste\] using full General Relativity. Independently of the present results and former work, we are able to prove that this strong growth is impossible if Bekenstein’s limit \[bek\] holds (P.S. Custódio and J.E. Horvath, in preparation).

## V Conclusions

We have analyzed in this work some features of the evolution of PBH mass in the radiation-dominated era. We found that within FRW models with null curvature there is no room for an explosive growth of their masses, which would need relativistic speeds of the gravitational radius. This result is in line with the previous analysis of Refs.\[ste\] and \[NOS\]. Indeed the very modest value of the maximum fractional gain ($`0.04`$) argues against the possibility discussed in Ref.\[Siv\]. This means that essentially all the methods applied to limit their abundance make sense since evidence for these objects could lead to very valuable insights. We have also addressed in some detail two ways to preclude evaporation: a large proper motion altering the absorption-emission properties of the semiclassical PBHs and the possibility that they could be ”dragged” and stay on the critical mass curve. In the first case we have given an ultrarrelativistic motion of the holes at their very birth ($`\gamma 36`$) ; and demonstrated that the second situation is one of unstable equilibrium, and therefore that PBHs do cross the critical mass curve as discussed in \[NOS\] to evaporate (Fig.1). Other forms to avoid the onset of the evaporation might in principle be possible, but are not known to us.

## VI Acknowledgements

Both authors wish to thank the São Paulo State Agency FAPESP for financial support through grants and fellowships. J.E.H. has been partially supported by CNPq (Brazil).

## VII Appendix: a general argument against a linear growth

We shall now discuss further consequences of a linear growth of the mass as derived in Ref. \[Siv\]. If we assume that it is correct and substitute $`M(t)=M_i(t/t_i`$ back into $`\dot{M}(t)=B\varrho _{rad}(t)M^2(t)`$ and impose that $`\varrho _{rad}(t)=\varrho _{rad}(t_i)(t_i/t)^2`$, we obtain for the r.h.s.

$$\frac{dM}{dt}=B\varrho _{rad}(t_i)=constant;$$
(50)
for arbitrary times. Using again the same expression $`M(t)=M_i(t/t_i)`$ and rearranging yields

$$\frac{M_i}{t_i}=BM_i^2\varrho _{rad}(t_i);$$
(51)
which shows that a linear solution would require a definite value of the radiation density at the formation time. Since initial black hole masses and the environment radiation density must be independent, we conclude that a linear growth is unphysical.

We can generalize the above fine tuning constraint for arbitrary conditions provided reasonable physics still holds. Suppose that we do not know the behaviour of the radiation density at all. Then, it would be reasonable to guess for it some function of the form $`\varrho _{rad}(t)=\varrho _{rad}(t_i)F(t/t_i)`$.

If we substitute this guess into eq.(17), we obtain

$$(11/u)=M_iB\varrho _{rad}(t_i)_1^u𝑑u^{}F(u^{});$$
(52)
where $`u=(t/t_i)`$. A power-law ansatz $`F(u)=u^m`$ for some $`m`$ gives

$$(11/u)=\frac{1}{(m+1)}[BM_i\varrho _{rad}(t_i)](u^{(m+1)}1).$$
(53)
The solution for $`m=2`$ displays the same unphysical fine tuning as before, since that it requires $`BM_i\varrho _{rad}(t_i)=1`$.

On general grounds there is no physical solution for $`m`$ different from $`2`$ which does not require some fine tuning between the initial mass $`M_i`$ and the radiation density at this moment. Therefore, we conclude that the solution $`M(t)t`$ is ruled out.

## VIII Figure captions

Figure 1. Evolution of PBHs in the absorption and evaporation regimes

The type of paths of a generic PBH are sketched in the PBH mass vs. cosmological time plane. The PBH is formed at the point marked with the star and undergoes a very small mass gain as long as its mass is larger than the critical mass (solid line), thus following an almost perfect horizontal path (long dashed) until the crossing time at the point marked with the cross. Since, as discussed in Section 3, the equilibrium is instantaneous, the PBH must drift to the evaporation region where the Hawking radiation will cause a decrease of the PBH mass and thus an evolution along the downwards bending path as indicated. Note that the mass as a function of time has a maximum on the critical mass curve; given that $`\dot{\sigma }`$ has the inverse sign than $`\ddot{M}`$. The path claimed in Ref. 3 is depicted with a dotted line and would need unlikely conditions as discussed in the text.


text:
# First-principles study of nucleation, growth, and interface structure of Fe/GaAs

## I Introduction

Magnetic materials are widely anticipated to be integrated into semiconductor-based microelectronics during the next decade or two.  A major component of this effort has focused on using ferromagnetic thin films as a source for creating spin-polarized electrical current in a semiconductor substrate, a process referred to as “spin injection.” Three criteria have been identified as important for technologically useful spin injection: (1) Substantial spin polarization of the injected current; (2) Electrons, rather than holes, serving as the spin-polarized carriers; (3) Curie temperatures for the source of order room temperature or higher.

Although several materials appear to meet one or two of these criteria, none has yet met all three. Several diluted magnetic semiconductors based on Be- and Mn-doping of ZnSe, and Mn-doping of GaAs, have recently demonstrated spin-injection efficiencies of greater than 50%, but only near liquid helium temperatures. Other magnetic semiconductors, including CdCrSe, can be doped both $`p`$\- and $`n`$-type and have Curie temperatures of order 100–200 K, but their performance as spin injection sources has yet to be examined.

In addition to the magnetic semiconductors, much current research continues to focus on one of the earliest studied candidate materials, Fe.  Besides offering the possibility of room-temperature injection of electron spins, Fe has the potential advantage of forming a nearly lattice-matched epitaxial film on an important semiconductor substrate, GaAs. However, despite this early promise and after considerable research investment, measured spin-injection efficiencies for Fe/GaAs remain frustratingly low, typically no larger than 1%.  The reason for this low efficiency is not definitively known, and much controversy surrounds its origin. A crucial issue, not yet settled, is whether the measured efficiencies reflect an intrinsic upper limit, or are simply due to technical limitations that may be overcome or circumvented. For example, much early effort focused on the nature of interface layers inferred to be “magnetically dead”, in the sense that they do not to contribute to the total magnetic moment of the Fe film. The appearance of dead layers is consistent with the formation (thermodynamically favorable in the bulk) of nonmagnetic FeAs complexes at an interface with the As-rich GaAs(001)-(2$`\times `$4) surface. This reasoning led to the development of two strategies for suppressing As diffusion into the Fe film: passivation of the As-rich surface by a surfactant such as sulfur; and growing on the Ga-rich GaAs(001)-(4$`\times `$6) surface (or, equivalently, As decapping prior to growth). Both approaches lead to films with magnetization onset in the range 4–8 monolayers and with essentially the full moment per Fe atom in all layers.

Notwithstanding the successful elimination of dead layers, measured spin-injection efficiencies have, until recently, remained at or below the 1% level. Recently, Schmidt et al. have argued that a more fundamental limitation exists for spin injection from a ferromagnetic metal into a semiconductor. They have shown that in the purely diffusive regime (where spins are scattered much less frequently than electrons) the spin-injection efficiency from a ferromagnet (fm) into a semiconductor (sc) is proportional to the ratio of their conductivities, $`\sigma _{\mathrm{sc}}/\sigma _{\mathrm{fm}}`$. For ferromagnetic metals this ratio is of order 10<sup>-4</sup>, and for typical device geometries suggests maximum injection efficiencies of 1% or less.

There remain several possibilities for circumventing this limitation on Fe sources. The first is to operate in the ballistic regime, where the contact resistance due to elastic backscattering at the interface (Sharvin resistance) will generally reduce the metal-semiconductor conductivity mismatch. Tang et al. have used concepts from mesoscopic transport to model injection into a two-dimensional electron gas, and find clear evidence for ballistic spin transport that otherwise vanishes in the diffusive limit.

A second possibility, relevant to the diffusive regime, is the use of tunnel contacts at the Fe/GaAs interface, which are expected to substantially reduce the conductivity mismatch and thereby increase the injection efficiency. Related research avenues concern the role of intrinsic Schottky barriers in controlling the spin-dependent tunneling through an interface. Indeed, Zhu et al. have recently demonstrated injection from Fe into GaAs with an efficiency of about 2% at room temperature; they suggest that the Schottky barrier formed between Fe and GaAs leads to a tunnel contact, thus circumventing the problems of conductance mismatch. Very recently, Hanbicki et al. have used a Schottky contact to an AlGaAs overlayer to inject spin into GaAs with efficiencies of 30% at low temperature and 9% at room temperature.

Since the prospects and limitations for spin injection from Fe into GaAs remain uncertain, we believe that further progress may benefit from a first-principles theoretical description of the interface. Our focus will be on the atomic structure of the interface and its resultant magnetic character, especially during the first few monolayers of growth. Theoretical studies of buried interfaces are notoriously difficult, for several reasons. First, experimental probes can provide only indirect information about physical and electronic structure, and hence are of limited utility for guiding theories. Moreover, real interfaces may—even after careful annealing—have atomic geometries very different from the ground state, hence the interface structure may depend on the precise growth history. In principle, one approach to this dilemma would be first-principles finite-temperature molecular dynamics simulation of film growth, which would properly account for the roles of deposition rate, surface diffusion, and incorporation into the substrate. In practice, molecular dynamics using density functional theory can at best simulate processes for $``$100 picoseconds—many orders of magnitude short of the experimentally relevant time scales, which may be milliseconds or longer. Kratzer and Scheffler have recently addressed this problem using a “first-principles kinetic Monte Carlo” method. However, the complexity of applying this method to the growth of Fe on GaAs—with three atomic species and many different possible processes—makes such an approach not yet feasible.

In this paper we approach the problem from two different limits. In Section II we consider the initial stages of interface formation, beginning with adsorption of isolated Fe adatoms on a bare GaAs substrate. We focus here on the potential-energy surface governing surface diffusion, and show that Fe-As chemistry may play a decisive role in the submonolayer regime. In Section III we consider a different limit, namely the thermodynamic ground-state interface structure of Fe films several monolayers thick. Here we concentrate on the magnetic character of the interface as a function of film thickness, and propose a mechanism that accounts for the experimentally observed delayed onset of ferromagnetic order for films of a few monolayers. In Section IV we discuss ways to reconcile our results from these two limiting regimes, and we propose future directions for further study.

## II Initial stages of growth

Our goal in this section is to identify and quantify those structural, magnetic, and chemical features that may be generic to the growth of Fe on GaAs. We do not attempt a definitive treatment of growth on a particular GaAs reconstruction or at a particular temperature. Instead, we focus on the following four questions: (1) Does the initial adsorption of atomic Fe prefer metallic bonding at highly coordinated sites, or does covalent bonding to Ga or As prevail? (2) Does adsorbed Fe nucleate the formation of compact islands, or do adsorbate-substrate interactions dominate the growth? (3) Is surface diffusion of Fe likely to be significant at growth temperatures? (3) Are the Fe magnetic moments in the submonolayer regime governed by strong Hund coupling (atomic moments) or by itinerancy effects (bulk moments)?

We begin by studying Fe adsorption on a fictitious, but chemically reasonable, surface of GaAs(001): a ($`2\times 1`$) reconstruction consisting of bulk GaAs terminated by a dimerized Ga layer. This fictitious surface is analogous to the dimerized $`\beta (2\times 4)`$ and $`\beta 2(2\times 4)`$ As-rich surfaces, but is quite different from the more complex $`\zeta (4\times 2)`$ surface believed to be the basis for the $`c(8\times 2)`$ reconstruction observed under Ga-rich conditions.

We use density-functional theory together with the generalized gradient approximation (GGA) for the exchange-correlation functional. For the calculations in this section we used Troullier-Martins pseudopotentials and a plane-wave basis with a kinetic-energy cutoff of 50 Ry, as implemented in the fhi96md code. The sampling of the surface Brillouin zone was equivalent to using 64 $`k`$-points for a ($`1\times 1`$) surface unit cell.

To construct the potential-energy surface (PES) for adsorption of atomic Fe on the clean surface, we computed the total energy, $`E_t`$, as a function of the adsorbate position $`(X,Y)`$ within the surface unit cell. To minimize interactions of the adsorbate with its periodic images we used a ($`2\times 2`$) supercell. For each adsorbate position $`(X,Y)`$, we fully relaxed the $`Z`$ coordinate of the adsorbate and the positions of all Ga and As atoms in the top two layers of the surface. This procedure was repeated for adsorbate positions sufficient to sample the PES with a resolution of about 0.25 Å. For many adsorbate positions, geometries with subsurface adsorption or with atomic positions exchanged were considered as well; in each case, the lowest-energy configuration was used to define $`E_t(X,Y)`$.

In order to illuminate the detailed role—if any—of magnetism in the adsorption energetics, we begin by first calculating the PES without allowing for spin-polarization. We then recalculate the PES while allowing for spin-polarization, and examine the differences that arise. If the adsorbate moment were completely localized and did not interact with the substrate, these two energy surfaces would be identical except for an overall shift of the energies. Thus, the differences that we find (described below) directly reflect site-dependent magnetic interactions between the Fe moment and the GaAs substrate.

Fig. 1 shows the PES calculated without spin polarization, for a portion of the surface unit cell shown in Fig. 2. We restrict the plot to the vicinity of the dimer rows because the energy landscape between dimer rows is quite flat and considerably higher in energy. Along the dimer rows the PES is periodic with repeat length $`a/\sqrt{2}`$; Fig. 1 shows slightly more than one full period along this direction.

Each point on this PES corresponds to a different relaxed geometry, determined solely by the in-plane adsorbate position $`(X,Y)`$. We focus on the four points marked in Fig. 1, which correspond to the four geometries shown in Fig. 2. Point (a), in the center of the PES, corresponds to adsorption at the high-symmetry pedestal site, midway between two dimers along a dimer row. Although this site is highly coordinated, leading to four Fe-Ga and two Fe-As bonds of nearly equal length (2.5 Å), it is energetically unfavorable. Indeed, this site is not even locally metastable: energy is gained by moving the adsorbate in any direction away from point (a).

The global minimum of the PES is a point (d) near the corner of the region shown in Fig. 1. This is 0.8 eV below point (a) and corresponds to the configuration of Fig. 2(d), in which the Fe adsorbate has partially “kicked out” one Ga atom from a surface dimer and taken its place. The large energy gain from this process suggests that the formation of such Fe-Ga heterodimers may act as a strong local trap for Fe, strongly suppressing surface diffusion of the Fe adatoms. If the kicked-out Ga atom subsequently diffuses away to a more stable adsorption site, this trapping effect will be further enhanced. Such effects, if not kinetically hindered, will bias the growth toward nucleating many small Fe islands.

The issue of kinetic barriers to forming Fe-Ga heterodimers can be analyzed by considering the possible routes leading to point (d) along the PES. For example, the reaction pathway from point (a) to (d) proceeds via point (c), which corresponds to the configuration shown in Fig. 2(c). This is the transition state (the highest energy configuration along the minimum-energy pathway) for the reaction taking (a) into (d), and is only 0.25 eV higher than the valley floor near (a). Part of this barrier arises from breaking the original Ga-Ga dimer bond, but apparently most of this energy cost is recovered by forming a more stable Fe-Ga heterodimer bond. Another 1.0 eV is recovered at point (d) by forming, in addition, a new Fe-As backbond to the substrate. The kicked-out Ga atom may remain bonded to the Fe atom, but it is energetically much more favorable to break this Fe-Ga bond (at a cost of 0.6 eV) and diffuse to a more favorable binding site—for a net gain up to 2.0 eV per Ga atom for the limiting case of incorporation into bulk Ga. Comparing these energy changes upon bond breaking, we conclude that Fe-As bonds are considerably more stable than Fe-Ga bonds.

The final adsorbate site we discuss is point (b), halfway between the pedestal site and the dimer bridge site. At this point, the Ga-Ga dimer bond has already broken to allow for more favorable Fe-Ga bonds to begin forming. Without spin polarization, the energy of this configuration is a local maximum, and represents a strong barrier to breaking Ga-Ga dimer bonds head-on. A similar barrier is found for a Fe adatom artificially constrained to approach the Ga-Ga dimer from directly above, as shown in in Fig. 3. For this adsorption route from above, an energy cost of $``$0.4 eV must be paid before the Ga-Ga dimer bond finally breaks; once paid, stable Fe adsorption can occur at the bridge site.

We turn now to the differences in the PES that arise from including spin polarization. We have repeated the calculation of the entire PES using spin-polarized GGA, allowing for any changes in the relaxed geometries. In general, we find only negligible changes to the relaxed atomic positions, and so the geometries shown in Fig. 2 continue to correspond to points on the new PES, shown in Fig. 4. Comparing the spin-unpolarized and spin-polarized energy landscapes (Figs. 1 and 4), several features deserve comment. The global minimum is the same for both, corresponding the Fe-Ga heterodimer formation via Ga “kick-out.” With spin polarization, the pedestal site is again unstable against adsorbate motion in any direction.

Differences are also apparent. The spin-polarized energy surface appears more corrugated; this is due almost entirely to a strong reduction of energies near point (b) (a local maximum without spin polarization) and the nearby dimer bridge site (a saddle point without spin polarization, but part of a low-energy trench with it). Thus by including spin-polarization, the barrier to head-on Ga-Ga dimer breaking is reduced to zero. The same difference is observed in dimer breaking by a Fe atom artificially constrained above the dimer: the barrier in Fig. 3 is eliminated by including spin polarization—demonstrating that this effect is not limited to surface diffusion.

To understand why allowing for spin-polarization changes some parts of the PES but not others, we consider two plausible explanations: chemical effects and magnetic effects. By the former we mean contributions to the total energy that depend primarily on the total valence electron density; by the latter we mean contributions related to the electron spin density. We focus on the configuration of Fig. 2(b), for which the change in the PES is particularly dramatic. In Fig. 5 we compare the spin-unpolarized and spin-polarized total valence-electron density in a plane containing the Fe adsorbate and its As and Ga neighbors. In both plots, Fe-As and Fe-Ga bonds are clearly visible, and the loss of both the Ga-Ga dimer bond and the Ga-As backbond is obvious. Most importantly, there are no large differences between the spin-polarized and unpolarized electron densities in this plane. We conclude that the changes in the PES due to spin-polarization cannot be attributed to changes in the valence electron density related to chemical bonds.

To evaluate the role of magnetic effects in determining the shape of the PES, we first consider how the total magnetic moment per cell varies with the Fe adsorbate position. In Fig. 6 we show the magnetic energy (the difference between spin-polarized and unpolarized energies), as a function of total magnetic moment, for about 50 different adsorption sites on the PES. Most of the sites have magnetic moments between 2 and 3 $`\mu _B`$, giving rise to magnetic energies between 0.4 and 0.8 eV. About ten sites have considerable higher moments, between 3.2 and 3.4 $`\mu _B`$, and correspondingly higher magnetic energies, between 1.0 and 1.6 eV. Thus these sites, all located near point (b) on the PES, give rise to large differences in the spin-polarized and unpolarized energy landscapes. These differences can be attributed to the development of unusually large magnetic moments for adsorption sites in the vicinity of energy barriers in the spin-unpolarized PES; these sites generally involve highly strained or partially broken bonds. The same trend is observed for the constrained adsorption from above (Fig.3). The magnetic energy gain at the configuration corresponding to the barrier is 2.5 eV, nearly twice as large as at the equilibrium adsorbate height. This difference again arises from the difference in magnetic moments: 3.9 $`\mu _B`$ at the barrier configuration vs. 3.2 $`\mu _B`$ at equilibrium. Thus we conclude that magnetic effects play an important role in determining those parts of the PES corresponding to highly strained or partially broken bonds.

We end this section by turning briefly to another process important during growth: the adsorption of Fe in the presence of preadsorbed Fe. We consider 1/4 monolayer of Fe in its lowest-energy configuration, i.e. incorporated as Fe-Ga heterodimers with (2$`\times `$2) periodicity. To simplify the discussion, we assume a starting surface from which the kicked-out Ga atoms have detached and diffused away. We then recalculate the (spin-polarized) potential-energy surface for this preadsorbed surface, again fully accounting for geometrical relaxation.

A portion of this PES is shown in Fig. 7(a), along with the initial preadsorbed surface structure in Fig. 7(b). The area shown is roughly twice that of Fig. 4, and a general comparison of the corresponding portions reveals the following. First, the overall corrugation of the preadsorbed PES is smaller than for the clean surface, i.e. preadsorbed Fe lowers the barriers to surface diffusion. The locations on this PES labeled (e) through (h) correspond to global minima on the PES of the clean surface. Point (e), located at the position of the preadsorbed Fe, is here a local energy maximum. This suggests that Fe-Fe bonding is not yet favorable in this low coverage regime, at least relative to the further formation of Fe-Ga heterodimers. Points (f) through (h) are local minima, indicating that formation of Fe-Ga heterodimers remains favorable for the preadsorbed surface. The global minimum, point (h), corresponds to the staggered arrangement of Fe-Ga heterodimers shown in Fig. 7(c). Occupying all of these adsorption sites results in an Fe coverage of 1/2 monolayer.

Although we do not explicitly calculate energy surfaces for higher coverages, we can make plausible inferences for their basic features based on the PES of Fig. 7. Even for a surface completely terminated by Fe-Ga heterodimers, the strong preference for Fe to form backbonds to As—as evidenced by the local minimum labeled (f) in Fig. 7(a)—should persist. Thus we speculate that continued deposition of Fe will ultimately lead to Fe-Fe dimers atop the As-terminated substrate. An obvious consequence of such a configuration would be the release of Ga (a full monolayer in this surface model). Moreover, in Section III we will show that for Fe films of several monolayers, a floating Ga (or As) adlayer can act as a surfactant, lowering the surface energy by as much as 2 eV/(1$`\times `$1) cell.

To conclude this section, we have shown that the initial stages of adsorption of Fe are dominated by strong local Fe-As chemistry. For a surface terminated by Ga, this chemistry leads to facile breaking of surface Ga dimers by sequential kicking-out of the two Ga atoms, first forming Fe-Ga heterodimers and, finally, surface Fe dimers bonded to subsurface As. Excess Ga may be released as a result of this rebonding mechanism. Magnetism plays an auxiliary role in the process by lowering the potential-energy barriers to breaking apart surface Ga dimers. Since the formation of Fe-As bonds leads to efficient trapping of Fe—especially if the excess Ga diffuses away—surface diffusion of Fe will probably be strongly suppressed.

## III Ground-state interface structure

As a complementary approach to the study of single Fe adsorption and diffusion, we have studied epitaxial interfaces of Fe/GaAs. In this section, our focus is threefold: (1) to determine the stability and magnetic character of Fe films for different interface structures and film thickness; (2) to investigate the role of As or Ga adlayers on Fe films, in particular the extent to which they may account for the observed out-diffusion of substrate atoms to the surface; (3) to examine the possibility of antiferromagnetic order as the origin of the observed magnetic quenching of Fe films with thicknesses of just a few monolayers.

The Fe/GaAs interface structures considered in this section are shown in Fig. 8 for the the As-terminated GaAs(001) interface. Model A is an atomically abrupt interface of bcc Fe and zinc-blende GaAs. Because the lattice constant of bcc Fe ($`a=2.866`$ Å) is almost half that of the substrate ($`a=5.654`$ Å), the epitaxial relationship is Fe(001)$`100`$ $`||`$ GaAs(001)$`100`$ and the atomic density of each Fe layer is twice that of the substrate layer, with a strain of only 1.3%. With respect to the bcc Fe lattice, there are vacancy sites in the adjacent GaAs lattice. Models B and C are built by filling these vacancy sites with Fe atoms one by one; we will refer to these interfaces as “partially intermixed” and “fully intermixed,” respectively. For each of these models, we considered film thicknesses ranging from 0.5 ML \[one Fe atom per $`(1\times 1)`$ of GaAs(001)\] up to 3.5 ML. We also considered three analogous models for the Ga-terminated interface. Both As- and Ga-terminated interface models may have relevance for experiments with the more commonly used Ga-rich surfaces: for example, we have already seen in the previous section that substitutional displacement of Ga on a Ga-rich surface leads to an As-terminated Fe/GaAs interface.

For the calculations in this section we again used the generalized gradient approximation (GGA), here with an ultrasoft pseudopotential for Fe. The electronic wave functions and densities were described by a planewave basis with cutoff energies of 16 Ry and 160 Ry, respectively. For each structural model, six or seven atomic layers were used for the substrate; the bottom As or Ga layer was passivated by pseudohydrogen atoms and fixed during the structural relaxation. The relative stability of each slab is given by its formation energy,

$$E_{\mathrm{form}}^{\mathrm{slab}}=E_t\underset{i}{}N_i\mu _i,$$
(1)
where $`E_t`$ is the total energy of the slab, $`N_i`$ the number of atoms of each chemical type, and $`\mu _i`$ their chemical potentials. To eliminate the contribution of the pseudohydrogen layer to the formation energy, we define the formation energy of each structural model with respect to a common reference structure,

$$E_{\mathrm{form}}^{\mathrm{model}}=E_{\mathrm{form}}^{\mathrm{slab}}E_{\mathrm{form}}^{\mathrm{ref}}/2.$$
(2)
The reference structure is an ideal GaAs slab passivated on both sides (hence the factor 1/2) by pseudohydrogen atoms. Assuming the thermodynamic equilibrium condition, $`\mu _{\mathrm{GaAs}}=\mu _{\mathrm{As}}+\mu _{\mathrm{Ga}}`$, and taking $`\mu _{\mathrm{GaAs}}`$ and $`\mu _{\mathrm{Fe}}`$ from bulk structures, the formation energy is expressed as a function of $`\mu _{\mathrm{As}}`$ only within the thermodynamically allowed range $`\mu _{\mathrm{GaAs}}\mu _{\mathrm{Ga}(\mathrm{bulk})}<\mu _{\mathrm{As}}<\mu _{\mathrm{As}(\mathrm{bulk})}`$. Here, the lower (upper) limit corresponds to the Ga-rich (As-rich) environment, and $`\mu _{\mathrm{Ga}(\mathrm{bulk})}`$ and $`\mu _{\mathrm{As}(\mathrm{bulk})}`$ are again determined from their bulk structures.

Figure 9 shows the calculated formation energies as a function of Fe film thickness. We first discuss Fe films with the As-terminated interface. At a Fe coverage of 0.5 ML model C is most stable, while models A and B are $``$1 eV higher. This energetic ordering persists up to 1 ML coverage, but changes above that point: the formation energy of model C increases by 1.8 eV from 1 ML to 2 ML, while that of the A and B become nearly independent of coverage. After the formation of 2 ML, the formation energy of all three models does not change much by adding more Fe layers. Thus, the low-energy interface structure turns out to be model A—the abrupt interface. For Fe films with the Ga-terminated interface, a similar trend applies but with a different outcome. At 0.5 ML, models B and C are equally stable, and at 1 ML model C is in fact the most stable. With additional Fe layers, model C becomes very unstable, similar to the case of the As-terminated interface. At higher coverages, however, model B—the partially mixed interface—is most stable.

These results can be understood as a competition between maximizing the coordination of Fe atoms and minimizing the concentration of atoms in the interfacial region. At low Fe coverage ($`1`$ ML), model C is energetically most favored, because Fe atoms can maximize their coordination with substrate atoms. At higher coverages, where the interface is well defined, this arrangement becomes unstable relative to the less intermixed interfaces. The reason for this crossover is the two extra Fe atoms per ($`1\times 1`$) surface cell (relative to model A). The excess electrons from the extra Fe atoms fill antibonding orbitals and thus weaken the interface bonding. This weakened bonding is also evident in the relaxed interface separations: for both As- and Ga-terminated interfaces, the position of the first Fe monolayer is $``$0.8 Å higher in model C and $``$0.4 Å higher in model B, compared to model A (see the relaxed structures in Fig. 8). Similar reasoning explains why model B is more stable than model A for the Ga-terminated interface, while the opposite holds for the As-terminated interface.

The present results emphasize that the low-energy atomic structure at low Fe coverage may not be extrapolated to high coverages. This should be true also in real growth situations, since Fe atoms at low coverages would take positions one or two layers deep so as to maintain maximal coordination. Further deposition of Fe leads to an partially mixed interface that ultimately becomes unstable. Therefore, it is likely that substantial rearrangement of the atomic structure occurs during the film growth, provided the temperature is sufficiently high. This may partly account for the experimental observation  that substrate atoms, especially As atoms, diffuse out to the surface during the Fe growth even for room-temperature deposition.

To follow this reasoning, we examined the effect of a Ga or As adlayer on top of the Fe film. In Fig. 10 we compare the formation energies of four selected interface structures with and without an adlayer. Both Ga and As adlayers stabilize the surface substantially, regardless of the interface structure and film thickness. This is not surprising, since adlayers increase the coordination of surface Fe atoms. The energy gain from As adlayers is particularly large, amounting to more than 1.5 eV/$`(1\times 1)`$ in most cases. Adlayers on 1 ML of Fe show strong dependence on the adatom site: the low-energy site is the position extended from the substrate GaAs lattice, suggesting a strong covalent bonding between the substrate and the adlayer, mediated by the intervening Fe layer.

The results for adlayer structures, together with the interface energetics, provide strong theoretical evidence that substitutional adsorption and/or atomic exchange is an essential process during the Fe growth—a finding which was anticipated by experiments and which we explicitly demonstrated in Section II. By substitutional adsorption processes, Fe atoms can maximize their coordination at every stage of the growth. These processes result in the segregation of substrate atoms to the surface, and simultaneously facilitate optimal interfacial atomic densities and thus relatively stable interface structures. Our calculations have shown that such processes are inevitable and that the energy gain from them is quite substantial.

We now turn to magnetic properties of the interface structures. Calculated spin moments (per Fe atom) of various model structures are given in Figs. 10 and 11. For bare Fe films without Ga or As adlayers, average Fe spin moments are much enhanced compared to the calculated bulk value for Fe of 2.33 $`\mu _B`$. Fig. 11 shows that for films of 3 ML this enhancement is still substantial, and relatively insensitive to the specific interface structure or substrate termination. For Fe films with adlayers, Fig. 10 shows that both As- and Ga-adlayers suppress total magnetic moments, by as much as 1 $`\mu _B`$ in several cases.

In order to examine local variations in the total magnetic moments, we calculated partial moments within a sphere centered on each atom, in analogy to the muffin-tin sphere in all-electron approaches. We used a sphere radius of 1.2 Å, which is slightly shorter than the minimum bond length between atoms. Table I shows the results for three selected structures. Local Fe moments are significantly enhanced ($``$3.0 $`\mu _B`$) at the surface layer. For the relatively thin Fe films considered in this work, the buried and interfacial Fe layers also have sizable enhancement. The adlayer suppresses spin moments not only of the top-layer Fe atoms but also slightly of the second-layer Fe atoms. Both the enhancement and suppression can be understood in terms of changes in the Fe $`d`$-band width due to changes in coordination and symmetry. On the other hand, substrate atoms have small induced moments up to 3–4 layers deep, always with negative spin moments for the interfacial atoms. The sign of spin moments changes to positive at deeper layers, suggesting the formation of a spin-density wave.

Recent magnetic measurements have observed a delayed onset of the magnetic phase at $``$2 ML for Fe films grown on the Ga-rich GaAs(001)-$`(2\times 6)`$ and $`(4\times 2)`$ surfaces.  Of the various surface structures we have considered here, one shows a similar quenching effect. The As-terminated interface with 1 ML of Fe (the first structure in Fig. 10) has a zero net moment when capped by an As adlayer. Since this structure also has the lowest formation energy for an Fe coverage of 1 ML, it might be anticipated that the delayed onset is due to the quenching of magnetism from the strong covalent bonding between the As atoms via Fe atoms. However, there is another possibility, namely the formation of antiferromagnetic (AF) order. We have considered the possibility of AF order for the simplest structures with 1 ML of Fe in the same plane, so that one Fe atom has spin up and the other down. The results, summarized in Fig. 12, indicate that films with an As adlayer become more stable by 0.05–0.09 eV/$`(1\times 1)`$ upon the formation of AF order, while bare films and films with a Ga adlayer become unstable by 0.2–0.5 eV/$`(1\times 1)`$. We note that the lowest energy structure with quenched magnetism is unstable toward AF order. We propose on this basis that the observed delayed magnetic onset is due to the initial formation of AF order, and tightly correlated with the out-diffusion of As atoms to the surface.

While the microscopic origins driving the formation of AF order are not clear in detail, it appears that the surface As atoms mediate the AF order between Fe atoms. A related observation is that the AF-stabilized surface structure is locally similar to the tetragonal Fe<sub>2</sub>As structure in the AF ground state: the unit cell consists of bimolecular units with two Fe atoms in the same plane, plus one As and one Fe atom on each side of the Fe plane.  Thus it is plausible that AF order in Fe/GaAs interfaces at low coverage arises from the formation of the Fe<sub>2</sub>As-like structures. This local similarity is broken by additional Fe adsorption, so that ferromagnetic order is ultimately favored.

To conclude this section, we have shown that the low-energy interface structure is different for low and high Fe coverages. The addition of Ga- or, especially, As-adlayers substantially lowers the surface formation, providing a theoretical basis for the experimentally observed out-diffusion of substrate atoms to the surface. We have shown that Fe films on the GaAs substrate usually assume a ferromagnetic ground state even at very low coverages, while for certain cases an As adlayer can induce antiferromagnetic order in the Fe film. Since the latter cases generally have lower formation energies and are thus likely to form, we propose that the observed magnetic quenching for very low Fe coverage may be due to the formation of antiferromagnetic order.

## IV Summary and future directions

We have taken two complementary approaches to understanding the nucleation and growth of Fe on GaAs—one focusing on the behavior of single Fe adatoms deposited on clean and partially preadsorbed GaAs(001), the other focusing on the interface structure of complete Fe films at coverages up to several monolayers. Although a detailed growth history of Fe/GaAs interfaces cannot yet be described, our studies suggest four generic features that may play an important role.

First, for low coverages we have identified a very strong driving force for Fe to be highly coordinated. For single adatoms, this tendency is strong enough to spontaneously break surface Ga-Ga and Ga-As bonds in order to form Fe-As bonds. Although magnetic effects do not dominate this chemistry, they play an interesting auxiliary role by essentially removing energetic barriers from the reaction pathways. For half-monolayer films constrained by symmetry, these local surface chemical reactions are not possible, and the Fe atoms respond by occupying subsurface sites with high coordination to either Ga or As.

A second generic feature is the crossover from a preference for strongly intermixed films to less intermixed or even abrupt films. This crossover occurs between one and two ML for both As- and Ga-terminated interfaces. Its origin is not magnetic: the films exhibit moments per atom larger than the bulk value even for very low coverages, and the moments converge to their bulk value long after the crossover occurs. Rather, it arises from the competition between maximizing the coordination of Fe atoms (which favors intermixing) and minimizing the amount of excess interfacial Fe (which favors abrupt interfaces). For well-defined Fe films of 2 ML or more, the latter effect dominates and sharper interfaces become energetically preferred.

A third finding, common to all interfaces we have studied, is that Ga- and As-adlayers dramatically reduce the formation energies of Fe films. This stabilizing effect occurs for both Ga- and As-terminated interfaces, for both intermixed and abrupt interfaces, and for all film thicknesses considered. It is especially striking for As-adlayers, which can reduce the film formation by as much as 50%. We also find that this stabilization is generally accompanied by a suppression of the total magnetic moment of the film; since this is due to reduced local moments in the topmost layer or two, the effect is largest for thinner films. Indeed, for 1-ML Fe films with an As-adlayer, antiferromagnetic order can be more stable than ferromagnetic order.

Taken together, these three generic features collectively imply a fourth: the diffusion of Ga or As atoms from the interface to the surface of the Fe film. We have shown explicitly how Ga can be released by adsorption of individual Fe adatoms and the subsequent “kick-out” of Ga from surface dimers. Even for GaAs surface reconstructions that do not consist of Ga dimers, we speculate that the same strong Fe-As chemistry would again lead to the release of surface Ga atoms. We have also shown that during the growth—between one and two ML—a spontaneous rearrangement of the interface morphology is likely to occur, again leading to the release of either Ga or As. Although we do not speculate about the details of this atomic rearrangement, we have also shown that ultimately the liberated Ga or As is likely to play the role of a floating surfactant layer.

Finally, we mention a possible avenue for further research. One difficulty with theoretical studies of interface is the paucity of macroscopic observables that can be directly related to the microscopic interface structure. Schottky barriers are extremely sensitive probes of interface structure, varying by as much 25% for local changes in interface geometry. Schottky barrier heights have been measured in Fe/GaAs interfaces to be of order 0.7 eV, and thus represent a useful probe of interface microstructure. Moreover, Schottky barriers may actually be a necessary ingredient for circumventing the intrinsic limitations on spin injection from a ferromagnetic metal into a semiconductor. A theoretical understanding of their dependence on Fe/GaAs interface structure—including substrate termination and reconstruction, degree of intermixing, and magnetic character—would be a great asset.

## V Acknowledgments

This work was funded in part by ONR and in part by the Deutsche Forschungsgemeinschaft. Computational work was supported in part by a grant of HPC time from the DoD Major Shared Resource Center ASCWP. S.C.E. and S.H.L. acknowledge generous support from the Alexander von Humboldt Foundation.