text:
# Incomplete information and fractal phase space

## Abstract

The incomplete statistics for complex systems is characterized by a so called incompleteness parameter $`\omega `$ which equals unity when information is completely accessible to our treatment. This paper is devoted to the discussion of the incompleteness of accessible information and of the physical signification of $`\omega `$ on the basis of fractal phase space. $`\omega `$ is shown to be proportional to the fractal dimension of the phase space and can be linked to the phase volume expansion and information growth during the scale refining process.

PACS number : 02.50.Cw,02.70.Rr,89.70.+c,89.75.Da

Introduction

A fundamental hypothesis in Boltzmann statistical physics, explicit or not, is that the probability distribution is complete so that $`_{i=1}^wp_i=1`$, where $`w`$ is the number of all the possible states labelled by $`i`$ and $`p_i`$ is the probability that the system is at the state $`i`$. This is the complete probability normalization. In other words, one can treat a complete or asymptotically complete set of states of the system under consideration. This hypothesis of complete statistics implies that the following two conditions must be satisfied. 1) All possible states of the system of interest are well known, which requires a complete knowledge of the dynamics of the system. 2) The known states are accountable and calculable, or in other words, the information is accessible to us. Are these two conditions always satisfied?

As the study of complexity advanced, it is revealed that chaotic and fractal behaviors are ubiquitous in nature, and that a fractal phase space in general cannot be completely and exactly exploited. Any calculation based only on the accessible states or the differentiable (integrable) points in this phase space is necessarily incomplete due to the rejected or unknown (singular or unaccessible) points. So the complete descriptions may be satisfactory if and only if these unaccessible points are negligible. If it is not the case, these points should be taken into account or the theoretical descriptions would be aberrant. The question is : how to include the states of a physical system in a theory when they are not accessible to us?

Recently, we have developed a method in this direction : a statistical mechanics based upon the notion of incomplete information which has been shown capable of giving $`coherently`$ the Tsallis $`q`$-exponential distributions within the nonextensive statistical mechanics (NSM) and, moreover, useful quantum distributions for some systems of correlated electrons. In the present paper, I am presenting some arguments to support the hypothesis of information incompleteness. The physical significations of the incompleteness parameter will be discussed on the basis of fractal phase space.

Incomplete information hypothesis

In the conventional probabilistic science, it is supposed that our ignorance, or information, is $`completely`$ accessible and may be calculated by, e.g., the Hartley formula $`I_i=\mathrm{ln}(1/p_i)`$ and Shannon information measure $`I=_i^wp_i\mathrm{ln}(1/p_i)`$, as well as by other generalized information measures, under the harsh condition that all possible states are accessible to us so that all probabilities sum to one.

However, the hypothesis of incomplete information admits simply that a part of our ignorance about complex system may not be accessible to our treatment. It cannot be calculated from probability distributions. In other words, the information calculated by the above mentioned method on the basis of the complete probability normalization may be incomplete because here we only have $`_{i=1}^vp_i=\mathrm{\Omega }`$ where $`v`$ is the number of the accessible or accountable states and may be greater or smaller than $`w`$, the total number of states. $`\mathrm{\Omega }`$ represents the incompleteness of the treatment and necessarily linked to the nature of the system. Logically, the origine of this incompleteness may be attributed either to the partial knowledge of the dynamics or to the unaccessible (incalculable) states of the system.

Incomplete information hypothesis has been motivated originally by some fundamental problems encountered in NSM. It replaces the complete probability normalization by

$$\underset{i=1}{\overset{v}{}}p_i^\omega =1,$$
(1)
where $`\omega `$ is referred to as incompleteness parameter which equals unity if the probability distribution is complete. Then on the basis of the Hartley additive information measure $`I_i=\mathrm{ln}(1/p_i)`$ and the nonadditive generalization of Hartley formula $`(I_i)_\omega =\mathrm{ln}_\omega (1/p_i)=\frac{(1/p_i)^{\omega 1}1}{\omega 1}`$, two incomplete entropies $`S_{\omega 1}=k_ip_i^\omega \mathrm{ln}p_i`$ and $`S_{\omega 2}=k\frac{1_ip_i}{1\omega }`$ $`(\omega R)`$ are deduced. These formalisms not only give in a coherent way the nonextensive incomplete distributions for canonical ensemble $`p_i[1(1\omega )x]^{\omega /(1\omega )}`$ which has been proved very useful for many systems having non-gaussian distributions, they also lead to quantum distributions of which the first results of the applications to correlated electron systems seem quite interesting and promising.

With the incomplete normalization Eq.(1), the parameter $`\omega `$ can be uniquely related to the incompleteness $`\mathrm{\Omega }`$ by

$$\underset{i=1}{\overset{v1}{}}p_i^\omega +(\mathrm{\Omega }\underset{i=1}{\overset{v1}{}}p_i)^\omega =1,$$
(2)
for any $`0<p_i<1`$, $`v`$ and $`\mathrm{\Omega }`$. Figure 1 shows the relation between $`\omega `$ and $`\mathrm{\Omega }`$ for a given incomplete distribution $`\{p_{i=1\mathrm{}5}\}=\{0.1,0.15,0.2,0.25,\mathrm{\Omega }0.7\}`$. We see that, in general, $`\mathrm{\Omega }>1`$ and $`\mathrm{\Omega }<1`$ lead to $`\omega >1`$ and $`\omega <1`$, respectively. If $`\mathrm{\Omega }`$ is such that a certain $`p_k1`$ among $`\{p_{i=1\mathrm{}5}\}`$, then $`\omega \mathrm{}`$. If $`\mathrm{\Omega }0`$, $`\omega 0`$.

Incomplete normalization in its general form

In general, we should write

$$\underset{i=1}{\overset{v}{}}F_\omega (p_i)=1$$
(3)
in order to normalize incomplete distributions. The function $`F_\omega `$ should depend on the nature of the system and become identity function whenever information is supposed complete ($`\mathrm{\Omega }=1`$). The arithmetic average of a quantity $`𝐱`$ should now be given by

$$\overline{x}=\underset{i=1}{\overset{v}{}}F_\omega (p_i)x_i.$$
(4)
where $`x_i`$ is the value of $`𝐱`$ at the state $`i`$.

$`F_\omega `$ may be determined if the information measure and the distribution law are given. For example, with Hartley information measure and exponential distribution, $`F_\omega `$ can be showed to be identity function. In general, by entropy maximization through the functional

$$\delta [\underset{i}{}F_\omega (p_i)I(p_i)+\alpha \underset{i}{}F_\omega (p_i)+\beta \underset{i}{}F_\omega (p_i)x_i]=0$$
(5)
we get :

$$\frac{\mathrm{ln}F_\omega (p_i)}{p_i}=\frac{I/p_i}{I+\alpha +f_\omega ^1(p_i)}$$
(6)
or $`F_\omega (p_i)=C\mathrm{exp}[\frac{I/p_i}{I+\beta f_\omega ^1(p_i)}𝑑p_i]`$ where $`\alpha `$ and $`\beta `$ are the multipliers of Lagrange in Eq.(5), $`I(p_i)`$ is the information measure, $`p_i=f_\omega (x_i)`$ the distribution function depending on the parameter $`\omega `$, and $`C`$ the normalization constant of $`F_\omega `$.

Within this frame, we would be able to group many entropies into the family of incomplete entropies if we suppose that the Hartley formula $`I_i=\mathrm{ln}(1/p_i)`$ holds. For example, the additive entropy proposed by Gorban and coworkers to generalize Gibbs-Shannon entropy and to obtain the long and short tail distributions may be considered as the expectation of Hartley information with $`F_\alpha (p_i)=p_i\alpha (p_i+1)`$ for an incomplete description discussed by the authors. In this case, the incompleteness of the description may be given by $`_{i=1}^vp_i=\frac{1\alpha v}{1\alpha }`$. When $`\alpha =0`$, the description becomes complete and Gorban entropy recovers Gibbs-Shannon one. The second example is the entropy due to Burg for solving geophysical problems : $`S_i\mathrm{ln}p_i`$. To have this, we may put $`F_\omega (p_i)=1/v`$ which is a pseudo-equiprobability. Other examples are the entropies of Belis-Guiasu with $`F_\omega (p_i)=w_ip_i`$ for the study of the information in cybernetic systems and the entropies of Aczel-Daroczy with $`F_\omega (p_i)=p_i^\gamma `$. So the incomplete normalization seemingly provides a plausible approach to link certain special entropies to elementary information.

About incomplete information in fractal phase space

Through a detailed analysis of the entropy and information on fractal supports, Naschie shows a connection between the information and the topological dimension of some fractal sets. It is also shown that the dynamics of this kind of phase space is quasi-ergodic for $`n4`$ where $`n`$ is the dimension of the set or the number of particles in the fractal phase space of topological dimension $`d_f`$. The results of the analysis tell us that the Boltzmann entropy and Gibbs-Shannon entropy for Cantor set increase with increasing capacity dimension given by $`d_c^{(n)}=(1/d_f)^{n1}`$ which in turn increases with increasing dimension $`n`$ of the fractal set. On the other hand, $`d_c^{(n)}`$ should increase with decreasing topological dimension $`d_f`$. This behavior is similar to that of $`S_{\omega 1}`$ and $`S_{\omega 2}`$ with respect to the parameter $`\omega `$, which means that $`\omega `$ may have something to do with the topological dimension in fractal phase space.

In what follows, the information on fractal support will be studied from a different angle with the help of the incomplete normalization Eq.(1). First of all, we will show that this kind of power normalization is inevitable for the calculation of fractal information.

An incomplete normalization

For the sake of simplicity, let us consider a phase space in which the trajectory of a complex system forms a self-similar fractal structure, say, Sierpinski carpet (Figure 1). This means that the state point of the system can be found only on the black rectangular segments whose number is $`W_k=8^k`$ at $`k^{th}`$ iteration. Hence the total surface at this stage is given by $`S_k=W_ks_k`$ where $`s_k=l_0/3^k`$ is the surface of the segments at $`k^{th}`$ iteration and $`l_0`$ the length of side of the square space at $`0^{th}`$ iteration. If the segments do not have same surface, we should write $`S_k=_{i=1}^{W_k}s_k(i)`$.

Now let us suppose that the density of state is identical everywhere on the segments and that the dynamics in this phase space is quasi-ergodic : every point on the segments is equally visited, so that the probability for the system to be in the $`i^{th}`$ segment may be defined as usual by $`p_i=s_k(i)/S_k`$. This probability is obviously normalized. The problem is that, as discussed in , $`S_k`$ is an indefinite quantity as $`k\mathrm{}`$ or as the time $`t\mathrm{}`$. So strictly speaking, it can not be used to define exact probability definition or all the probabilities would be null or lose its additivity and normalization according to the countable additivity axiom of probability or total probability axiom which asserts that probabilities can be additive on a set of finite cardinals (number of elements). In addition, $`S_k`$ is not differentiable and contains inaccessible points. Thus the probability defined above makes no sense.

Alternatively, the probability may be reasonably defined on a integrable and differentiable support, say, the Euclidean space containing the fractal structure. To see how to do this, we write $`S_k=l_0^2(\frac{1}{3^k})^{dd_f}`$ for identical segments or, for segments of variable size,

$$\underset{i=1}{\overset{W_k}{}}[\frac{s_k(i)}{S_0}]^{d_f/d}=1$$
(7)
where $`S_0=l_0^d`$ (here $`d=2`$ for Sierpinski carpet) is a characteristic volume of the fractal structure embedded in a $`d`$-dimension Euclidean space, $`d_f=\frac{\mathrm{ln}n}{\mathrm{ln}m}`$ is the fractal dimension. $`n=8`$ is the number of segments replacing a segment of the precedent iteration and $`m=3`$ the scale factor of the iterations.

The microcanonical probability distribution at the $`k^{th}`$ iteration can be defined as

$$p_i=\frac{s_k(i)}{S_0}$$

so that

$$\underset{i=1}{\overset{W_k}{}}p_i^{d_f/d}=1.$$

This is just Eq.(1) with $`\omega =d_f/d`$. The complete probability normalization $`_{i=1}^{W_k}p_i=1`$ can be recovered when $`d_f=d`$.

It should be noticed that, in Eq.(7), the sum over all the $`W_k`$ segments at the $`k^{th}`$ iteration does not mean the sum over all possible states of the system under consideration. This is because that the segment surface $`s_k(i)`$ does not represent the real number of state points on the segment which, as expected for any self-similar structure, evolves with $`k`$ just as $`S_k`$. So at any given order $`k`$, the complete summation over all possible segments is not a complete summation over all possible states. But in any case, whatever is $`k`$, Eq.(7) and $`_{i=1}^{W_k}p_i^\omega =1`$ always hold for $`\omega =d_f/d`$.

In this simple case with self-similar fractal structure, the incompleteness of the normalization Eq.(1) is measured by the parameter $`\omega =d_f/d`$. $`d_f>d`$ means that the system has more states than $`W_k`$, the number of accessible states at given $`k`$. If $`d_f<d`$, the number of states is less than $`W_k`$. When $`d_f=d`$, the summation is complete at any order $`k`$, corresponding to complete information calculation.

$`\omega `$ and phase space expansion

Now we will discuss in a detailed way the incompleteness parameter $`\omega `$ and its physical meanings. Let us discuss this still in the case of self-similar fractal phase space.

As discussed in the case of chaotic phase space, $`\omega =\mathrm{ln}n/d\mathrm{ln}m`$ gives a measure of the incompleteness of the state counting in the $`d`$-dimension phase space. $`\omega =1`$ means $`d_f=d`$ or $`n=m^d`$. In other word, at the $`k^{th}`$ iteration, a segment of volume $`s_k`$ is completely covered (replaced) by $`n`$ segments of volume $`s_{k+1}=s_k/m^d`$. So the summation over all segments is equivalent to the sum over all possible states, making it possible to calculate complete information.

When $`\omega >1`$ (or $`\omega <1`$), $`n>m^d`$ (or $`n<m^d`$) and $`s_k`$ is replaced by $`n`$ segments whose total volume is more (or less) than $`s_k`$. So there is expansion (or negative expansion) of state volume when we refine the phase space scale. An estimation of this expansion at each scale refinement can be given by the ratio $`r=\frac{ns_{k+1}s_k}{s_k}=\frac{n}{m^d}1=(\frac{1}{m^d})^{1\omega }1=(\omega 1)\frac{(m^d)^{\omega 1}1}{\omega 1}`$. $`r`$ describes how much accessible states increase at each step of the iteration or of the refinement of phase space. The physical content of $`\omega `$ is clear if we note that $`\omega >1`$ and $`\omega <1`$ correspond to an expansion ($`r>0`$) and a negative expansion ($`r<0`$), respectively, of the the accessible state volume at each step of the iteration.

When $`\omega =0`$, we have $`d_f=0`$ and $`n=1`$, leading to $`r=\frac{1}{m^d}1`$. The iterate condition $`n1`$ means $`\omega 0`$, as proposed in references . $`\omega <0`$ is impossible since it means $`d_f<0`$ or $`n<1`$ which obviously makes no sense. We can also write : $`\omega 1=\mathrm{ln}(r+1)/\mathrm{ln}(m^d)=\mathrm{ln}(ns_{k+1}/s_k)/\mathrm{ln}(m^d)`$, which implies that it is the difference $`\omega 1`$ which directly measures the accessible state space expansion through the scale refinement.

$`\omega `$ and information growth

The expansion of the accessible state volume of a system in its phase space during the scale refinement should be interpreted as follows : the extra state points $`\mathrm{\Delta }=ns_{k+1}s_k`$ acquired at $`(k+1)^{th}`$ order iterate are just the number of unaccessible states at $`k^{th}`$ order with respect to $`(k+1)^{th}`$ order. $`\mathrm{\Delta }>0`$ (or $`\mathrm{\Delta }<0`$) means that we have counted less (or more) states at $`k^{th}`$ order than we should have done. $`\mathrm{\Delta }`$ contains the accessible information gain (AIG) through the $`(k+1)^{th}`$ iterate.

To illustrate the relation between this “hidden information” and the parameter $`\omega `$, let us suppose that the distribution is scale-invariant. At the iterate of order $`k`$, the average information contained on $`s_k`$ is given by $`I_k=_{s_k}p^\omega I(1/p)𝑑s`$. At $`k+1`$ order, $`I_{k+1}=_{ns_{k+1}}p^\omega I(1/p)𝑑s`$. Hence AIG is just $`\mathrm{\Delta }I=I_{k+1}I_k=_{(ns_{k+1}s_k)}p^\omega I(1/p)𝑑s=\sigma _I\mathrm{\Delta }`$, where $`\sigma _I=p^\omega I(1/p)`$ is the information density or the average information carried by each state. The relative AIG is given by $`\mathrm{\Delta }I/I_k=r=(1\omega )\frac{(1/m^d)^{1\omega }1}{1\omega }`$ which is independent of scale but dependent on scale changes. For given scaling factor $`m`$, the magnitude of $`\mathrm{\Delta }I`$ or $`r`$ increases with increasing difference $`|1\omega |`$. The sign of $`r`$ (or AIG) was discussed earlier. For given $`\omega `$, $`|\mathrm{\Delta }I|`$ increases with decreasing scaling. For $`\omega =1`$ or $`m=1`$, there is no information gain, corresponding to the case of complete information. According to the above discussions, the incompleteness parameter $`\omega `$ may be considered as a measure of the degree of chaos or fractal.

Conclusion

Summing up, we have argued that the information of a complex system we deal with may be incomplete because we only partially know the dynamics of the system or a part of the states of the systems is not accessible to our calculation. The parameterized normalization $`_ip_i^\omega =1`$ is justified on fractal support, where the incompleteness parameter $`\omega `$ is shown to be proportional to the fractal dimension of the phase space and can be directly linked to the phase volume expansion and information growth during the scale refining process.

Acknowledgments

I acknowledge with great pleasure the very useful discussions with Professor M.S. El Naschie.


format: text

text:
# From the bulk to monatomic wires: An ab-initio study of magnetism in Co systems with various dimensionality

## Abstract

A systematic ab-initio study within the framework of the local-spin-density approximation including spin-orbit coupling and an orbital-polarization term is performed for the spin and orbital moments and for the X-ray magnetic circular dichroism (XMCD) spectra in hcp Co, in a Pt supported and a free standing Co monolayer, and in a Pt supported and a free standing monatomic Co wire. When including the orbital-polarization term, the orbital moments increase drastically when going to lower dimensionality, and there is an increasing asymmetry between the $`L_2`$ and $`L_3`$ XMCD signal. It is shown that spin and orbital moments can be obtained with good accuracy from the XMCD spectra via the sum rules. The $`T_z`$ term of the spin sum rule is surprisingly small for the wires, and the reason for this is discussed.

The modern methods to prepare nanostructured systems made it possible to investigate the influence of dimensionality on the magnetic properties. A central question thereby is how the qualitative behavior will change when going from two-dimensional to one-dimensional systems because it has been predicted that there is no long-range magnetic order at finite temperature in infinitely extended one-dimensional systems with short-range magnetic interactions. In the past, several experimental investigations of monolayer nanostripes of Fe on vicinal surfaces of W Elmers et al. (1994); Hauschild et al. (1998) or Cu Shen et al. (1997) have been performed, with a stripe width down to 1-10 nm. Most recently, Gambardella et al. Gambardella et al. (2000) succeeded to prepare a high density of parallel atomic chains along steps by growing Co on a high purity Pt(997) vicinal surface in a narrow temperature range. The magnetism of the Co wires was investigated Gambardella et al. (2002) by the X-ray magnetic circular dichroism (XMCD). Below a blocking temperature a long-range magnetic ordering owing to the presence of anisotropy barriers was found on the time scale of the experiment. Applying a simple model of exchange coupled superparamagnetic clusters com , the anisotropy energy could be obtained from the shape of the magnetization curve above the blocking temperature and it appeared to be much larger than the one for a Co monolayer on Pt which — in turn — is much larger than the one of hcp Co. Accordingly, a large orbital moment of $`0.68\mu _\mathrm{B}`$ per Co atom was found, the highest value ever reported for a 3d itinerant electron system.

So far, magnetism in quasi-one-dimensional systems was studied mainly in the insulating material $`\mathrm{CsNiF}_3`$ where the magnetic Ni ions are arranged along linear chains which are well separated from each other so that the interchain interaction is only $`10^3`$ or less of the intrachain interaction. Because of the one-dimensional character of this spin system and the easy-plane anisotropy, magnetic solitons play an important role for the dynamical and thermodynamical behavior which has been investigated by neutron scattering experiments Mikeska and Steiner (1991). The discovery of magnetism in one-dimensional metallic systems opens up the chance to extend the research in many respects. First, by considering various vicinal surfaces the distance between the steps and hence the chains can be modified so that the transition from the well-isolated chains to interacting chains can be studied. Second, it is possible to grow, e.g., biatomic wires along the steps and to manipulate the wire length, and to study the respective influence on magnetic properties. Finally, we expect that the damping of magnetic excitations is larger for metals than for insulators and especially large for one-dimensional metals, possibly leading to peculiar properties of the nonlinear spin excitations Hołyst et al. (1986).

The already existing ab-initio calculations for monatomic transition metal wires considered magnetic moments Weinert and Freeman (1983); Bihlmayer et al. ; Dorantes-Dávila and Pastor (1998); Zhou et al. (1999); Spišák and Hafner (2002), exchange couplings Spišák and Hafner (2002) and magnetic anisotropies Dorantes-Dávila and Pastor (1998); Zhou et al. (1999); Eisenbach et al. (2002). It turned out Zhou et al. (1999) that in these one-dimensional systems the orbital correlations were essential and that the magnetic anisotropy and hence also the orbital moments could not be calculated reliably by the local-spin-density approximation (LSDA) in combination with spin-orbit (SO) coupling. Instead, a term correcting explicitly for the orbital correlations has to be added. It is one of the objectives of the present paper to work out the growing importance of the orbital correlations with decreasing dimensionality by calculating orbital (and spin) moments by the LSDA Perdew and Wang (1992) including SO coupling (in a self-consistent manner) without and with the orbital-polarization term Eriksson et al. (1990) which takes into account at least in part the orbital-correlation effects. The calculations are performed for Co atoms in hcp Co, in a Pt supported and in a free standing monolayer, and in a Pt supported and in a free standing monatomic wire. Another objective is to calculate the respective XMCD spectra and discuss the influence of dimensionality on the accuracy of the spin and orbital moments when these are obtained from the XMCD spectra via the sum rules Thole et al. (1992); Carra et al. (1993).

To focus on the pure effect of the dimensionality we fix the nearest-neighbor distance of the atoms for all considered systems to the one of fcc Pt (2.77Å). For the monolayers and wires we perform supercell calculations, i.e., large unit cells which model the structure are repeated periodically. In the case of the monolayer the supercell consists of two Pt and one Co {111} layers in the fcc stacking and a vacuum sheet corresponding to two empty layers. The vicinal Pt(997) surface with the Co wire at the steps is modeled by the supercell shown in Fig. 1 with two additional vacuum layers on the top.

The supercell of the free standing wire is constructed by removing all Pt atoms from this supercell. Test calculations have shown that the results change only very slightly when going to larger supercells. For the monolayers and the wires we assume perpendicular magnetization, for hcp Co it is along the c axis.

We use the tight-binding linear-muffin-tin-orbital (LMTO) method in the atomic sphere approximation Andersen and Jepsen (1984) in which we have implemented the SO coupling and the tools for calculating the XMCD spectra, and the WIEN97 code Blaha et al. (1990) which adopts the full-potential linearized-augmented-plane-wave (FLAPW) method Wimmer et al. (1981) in which the SO coupling and the tools for the calculation of the magnetooptical effects and XMCD spectra Kuneš et al. (2001a, b) and the orbital-polarization term Rodriguez et al. (2001) have been implemented. The magnetic orbital moments $`m_l=\mu _\mathrm{B}l_z`$ and spin moments $`m_s=\mu _\mathrm{B}\sigma _z`$ are calculated directly from the wave functions as well as via the sum rules Thole et al. (1992); Carra et al. (1993) from the absorption coefficients $`\mu ^+(ϵ)`$, $`\mu ^{}(ϵ)`$ and $`\mu ^0(ϵ)`$ for light with right-circular, left-circular and z-axis polarization at the $`L_2`$ and $`L_3`$ edge according to

$`l_z`$ $`=`$ $`{\displaystyle \frac{2I_\text{m}N_\text{h}}{I_\text{t}}},`$ (1)
$`\sigma _z`$ $`=`$ $`{\displaystyle \frac{3I_\text{s}N_\text{h}}{I_\text{t}}}7T_z,`$ (2)
$`I_\text{m}`$ $`=`$ $`{\displaystyle \underset{E_\text{F}}{\overset{E_\text{c}}{}}}\left[(\mu _\text{c})_{L_3}+(\mu _\text{c})_{L_2}\right]𝑑ϵ,`$ (3)
$`I_\text{s}`$ $`=`$ $`{\displaystyle \underset{E_\text{F}}{\overset{E_\text{c}}{}}}\left[(\mu _\text{c})_{L_3}2(\mu _\text{c})_{L_2}\right]𝑑ϵ,`$ (4)
$`I_\text{t}`$ $`=`$ $`{\displaystyle \underset{E_\text{F}}{\overset{E_\text{c}}{}}}\left[(\mu _\text{t})_{L_3}+(\mu _\text{t})_{L_2}\right]𝑑ϵ,`$ (5)
with the XMCD signal $`\mu _\text{c}=\mu ^+\mu ^{}`$ and with $`\mu _\text{t}=\mu ^++\mu ^{}+\mu ^0`$. Here $`N_\text{h}`$ is the number of d holes and $`T_z`$ is the expectation value of the magnetic dipolar operator

$$\widehat{T}_z=\frac{1}{2}[𝝈3\widehat{𝐫}(\widehat{𝐫}𝝈)]_z,$$
(6)
where $`𝝈`$ denotes the vector of the Pauli matrices. The quantities $`E_\text{F}`$ and $`E_\text{c}`$ denote the Fermi energy and a cutoff energy. For details for such type of calculations see Ref. Ederer et al., 2002. The $`T_z`$ term is negligible for cubic surroundings but it is expected to become more and more important when reducing the dimensionality of the system. Since it is very difficult to measure $`T_z`$, this term is often neglected in the spin sum rule when analyzing the experimental data. One of the objectives of this paper is to asses this critically.

The results of the LSDA calculations with SO coupling are given in the upper part of Table 1.

It should be noted again that thereby we use the same nearest-neighbor distance between the atoms, i.e., the one of fcc Pt, for all structures. As explained in Ref. Ederer et al., 2002, the spin and orbital moments derived experimentally via the sum rules correspond essentially to the respective part of the valence wavefunctions with 3d angular momentum, if the “background” due to all additional contributions to the experimental absorption spectra is subtracted appropriately. We therefore compare in Table 1 the moments as calculated directly from the 3d spin and orbital densities with those obtained from the sum rules when taking into account for the absorption spectra and for the $`T_z`$ term also only that part of the valence wave functions which has 3d character. For comparison, we give also the values of directly calculated moments including all angular momentum contributions.

For all the structures there is only a slight difference between the directly calculated moments with only d contributions and with all contributions. The spin moments increase only moderately with decreasing dimensionality. For the wire on Pt our directly calculated FLAPW spin value of $`2.06\mu _\text{B}`$ agrees well with the one of Ref. Bihlmayer et al., . There is also a satisfactory agreement between the directly calculated spin moments and those obtained from the spin sum rule when including the $`T_z`$ term. The $`T_z`$ contribution is negligible for hcp Co. In contrast, for the free and the Pt supported monolayers, the values of $`\sigma _z`$ obtained from the sum rule when neglecting $`T_z`$ are between 20% and 30% smaller than the directly calculated values. Astonishingly enough, when going to even stronger reduced dimensionality, i.e., for the Pt supported wire and even more for the free wire, the $`T_z`$ term again appears to be of minor importance. To figure out the reason for this startling result, Fig. 2 shows the energy distributionWu et al. (1994) for $`T_z`$ for all the structures, i.e., the integral over the energy resolved contribution to $`T_z`$ up to an energy $`E`$ (for the integrals up to the Fermi level, $`E=E_\text{F}`$, $`T_z`$ equals the quantity in eq. (2) ).

It becomes obvious that for the free standing monolayer and for the free standing wire, $`T_z`$ depends drastically on $`E`$ resulting in sharp structures in the energy distribution due to the two- or one-dimensional character of the systems. Going to the Pt supported monolayer and wire, the energy dependence of $`T_z`$ is less pronounced and the sharp structures are smeared out, indicating that these systems are not really two- or one-dimensional due to the presence of the substrate. However, $`T_z`$ obtained from the integral up to $`E_\text{F}`$ is very small for the free standing Co wire but larger for the other systems (except for hcp Co). This clearly demonstrates that the expected rule of thumb — a larger $`T_z`$ for the smaller dimensionality — does not necessarily hold. We have performed the same calculation for a free standing Ni wire. The energy distribution of $`T_z`$ is very similar but the Fermi level is shifted and $`T_z`$ is consequently much larger than for the Co wire. This again demonstrates that the symmetry arguments alone are not able to estimate the size of $`T_z`$.

The orbital moments obtained by the LSDA calculation with SO coupling are very similar (about $`0.15\mu _\text{B}`$) for hcp Co and for the Pt supported Co monolayer and wire, while the values are slightly larger (about $`0.18\mu _\text{B}`$) for the free monolayer and considerably larger for the free wire. Thereby, there is a good agreement between the directly calculated orbital moments and those obtained from the calculated spectra via the sum rules. The rather small values for the Pt supported monolayer ($`0.13\mu _\text{B}`$) and the Pt supported wire ($`0.14\mu _\text{B}`$) are in conflict with the large corresponding experimental valuesGambardella et al. (2002) ($`0.31\mu _\text{B}`$ and $`0.68\mu _\text{B}`$). As discussed in Ref. Zhou et al., 1999, the reason is presumably that LSDA + SO coupling does not appropriately account for orbital correlations. We therefore have redone all the calculations by including in addition the orbital-polarization term Eriksson et al. (1990). Table 1, lower part, gives the so obtained results for the orbital moments (the spin moments are nearly unaffected by this additional term). It should be noted that the calculations for hcp Co have been performed for the nearest neighbor distance of fcc Pt which explains the large moment of $`0.4\mu _\text{B}`$. For the Co monolayer on Pt (Co wire on Pt) our value $`0.43\mu _\text{B}`$ ($`0.92\mu _\text{B}`$) is even larger than the experimental value Gambardella et al. (2002) $`0.31\mu _\text{B}`$ ($`0.68\mu _\text{B}`$) which may be in part due to the fact that our calculations did not take into account any structural relaxations at the surface. For the free standing wire, the orbital moment appears to be extraordinarily high ($`2.3\mu _\text{B}`$). Obviously the experimental trend — increase of the orbital moment with decreasing dimensionality — is well reproduced when taking into account the orbital-polarization term. The increasing orbital moments result in an increasing asymmetry between the calculated $`L_2`$ and $`L_3`$ XMCD spectra (Fig. 3), in qualitative agreement with the respective experimental spectra (Fig. 2 of Ref. Gambardella et al., 2002).

Both, for the LSDA calculations with and without the orbital-polarization term, the directly calculated orbital moments agree very well with those obtained from the sum rule. Therefore, for Co in structures of various dimensionality the spin and orbital moments probably may be obtained with high accuracy from the sum rules (if including the $`T_z`$ term in the spin sum rule).

To conclude, we have shown that the orbital correlations have to be taken into account explicitly in order to reproduce the experimentally observed trend to higher orbital moments and larger asymmetries between the L3 and L2 XMCD spectra in Co systems with reduced dimensionality, and we are about to do this also for a systematic study of the magnetic anisotropy. In addition, we calculate ab initio the further parameters appearing in a model Hamiltonian for the monatomic wires, i.e., exchange interactions from the adiabatic spin wave spectra Grotheer et al. (2001), and the Gilbert damping factor Gilbert (1955). The final objective is to come to a comprehensive description of the thermodynamic properties (see Ref. com, ) and of the damped nonlinear excitations (along the lines of Ref. Hołyst et al., 1986) of monatomic magnetic linear chains.

The work at Brookhaven is supported by U.S. Department of Energy under Contract No. DE-AC02-98CH10886.


text:
# Instability of the rhodium magnetic moment as origin of the metamagnetic phase transition in 𝛼-FeRh

## I Introduction

In 1938 FallotFallot (1938); Fallot and Horcart (1939) discovered that ordered bcc Fe50Rh50 undergoes a first order metamagnetic transition from an antiferromagnetic (AF) ground state to a ferromagnetic (FM) phase with increasing temperature. This transition occurs at $`T_\mathrm{m}320\mathrm{K}`$ and is accompanied by a volume increase of about $`1\%`$.Kouvel (1966); Zakharov et al. (1964) The Curie temperature $`T_\mathrm{C}`$ of the FM phase is of the order of $`670\mathrm{K}`$.Zakharov et al. (1964); Kouvel and Hartelius (1962) In contradiction to the first hypothesis of Fallot and Horcart,Fallot and Horcart (1939) X-ray diffraction measurements showed that the transition is isostructural.de Bergevin and Muldawer (1961) From Mössbauer and neutron diffraction measurements one knows that the FM phase has collinear magnetic moments of $`3.2\mu _\mathrm{B}`$ per Fe atom and $`0.9\mu _\mathrm{B}`$ per Rh atom.Shirane et al. (1963) At low temperatures an AF-II spin structure is found with Fe moments of $`3.3\mu _\mathrm{B}`$ and with vanishing Rh moments (see Figure 1).Shirane et al. (1964); Kunitomi et al. (1971) Application of hydrostatic pressure suppresses the FM phase, i. e. for a critical pressure of $`60\mathrm{kbar}`$ the system immediately transforms from the AF to the paramagnetic (PM) phase.

An early explanation for this behavior was based on the phenomenological exchange inversion modelKittel (1960) by Kittel (which originally was designed to explain metamagnetic transitions in other materials like Cr-modified Mn2Sb with a layered magnetic structure). In this model the exchange parameter varies linearly with the lattice constant and changes sign for a critical value $`a_\mathrm{c}`$. However, experimental findings like the rather large entropy changeKouvel (1966); Zakharov et al. (1964); Annaorazov et al. (1996); Lommel (1969); McKinnon et al. (1970) at $`T_\mathrm{m}`$, which is of the order of $`\mathrm{\Delta }S(T_\mathrm{m})12.519.7\mathrm{J}\mathrm{kg}^1\mathrm{K}^1`$, as well as elastic properties could not correctly be described.McKinnon et al. (1970) Tu et al.Tu et al. (1969) used a different approach by considering the large difference of the low-temperature specific heat constants, $`\gamma `$, of the AF and the FM phases: Measurements suggest that $`\gamma _{\mathrm{FM}}5962.5\mathrm{mJ}\mathrm{kg}^1\mathrm{K}^2`$ is about four to six times larger than $`\gamma _{\mathrm{AF}}10.516\mathrm{mJ}\mathrm{kg}^1\mathrm{K}^2`$.Tu et al. (1969); Ivarsson et al. (1971) Based on these observations, Tu et al. explained the transition by a change in entropy of the band electrons between the AF and the FM phases. An estimation of the free energy shows that these contributions have the right order of magnitude to explain the AF–FM transition, if one assumes that the electronic densities of states at the Fermi level do not vary considerably from low temperatures up to the transition. However, since the FM phase cannot be stabilized at low temperatures for stoichiometric FeRh, $`\gamma _{\mathrm{FM}}`$ has been measured in slightly more iron-rich alloys. With respect to the strong sensitivity of the magnetic phase diagram to small departures from stoichiometry,Vinokurova et al. (1988) one may doubt that the low temperature specific heat can be considered to be independent of concentration.Ponomarev (1972); Moruzzi and Marcus (1992a) A further drawback considering the explanation of Tu et al. arises from the fact that adding $`5\%`$ Iridium boosts $`\gamma _{\mathrm{AF}}`$ by almost an order of magnitude to a value of $`\gamma _{\mathrm{AF}}=101\mathrm{mJ}\mathrm{kg}^1\mathrm{K}^2`$:Ivarsson et al. (1971); Fogarassy et al. (1972) The compound Fe49.5Rh45.5Ir5 also undergoes a metamagnetic transition with $`T_\mathrm{m}`$ shifted to higher temperatures.Kouvel (1966) Since the relation between $`\gamma _{\mathrm{AF}}`$ and $`\gamma _{\mathrm{FM}}`$ is reversed in this material, the previously sketched explanation cannot be applied to this case.

In 1992 Moruzzi and Marcus performed ab initio calculations using spin polarized density functional theory (DFT) in the framework of the local density approximation (LDA) and the augmented spherical wave (ASW) formalism.Moruzzi and Marcus (1992a); Moruzzi and Marcus (1992b); Moruzzi and Marcus (1993) They calculated the total energy of FeRh as a function of volume for different magnetic structures. They found that the AF-II spin structure is the ground state, whereas the FM structure is another stable solution with higher energy and a larger volume. For the moments they obtained $`\mu _{\mathrm{Fe}}=2.98\mu _\mathrm{B}`$, $`\mu _{\mathrm{Rh}}=1.02\mu _\mathrm{B}`$ (AF phase) and $`\mu _{\mathrm{Fe}}=3.15\mu _\mathrm{B}`$ (FM phase), in agreement with experimental results.

In this work, we present further ab initio total energy calculations of stoichiometric $`\alpha `$-FeRh using the ASW formalism and the generalized gradient approximation (GGA).Perdew et al. (1992) The corresponding code also allows to evaluate noncollinear alignment of spins (see Ref. Uhl et al., 1994 and references therein). For many cases LDA gives a reasonable description of the ground state properties of solids. However, in some cases it predicts the wrong ground state as e. g. for ironStixrude et al. (1994). For Fe, the hcp structure is found to be $`10`$mRy lower in energy than the bcc phase. Because FeRh has also a bcc ground state lattice structure (unlike pure Rh), it is necessary to go beyond LDA by using GGA. We have considered different possible magnetic ground state structures presented in section II. Based then on the specific energetic order of resulting energy versus volume curves, we discuss a new mechanism for the temperature driven metamagnetic transition in FeRh. In contrast to prior explanations, our model does not rely on ground state properties or low temperature data alone. Instead, we propose that thermal excitations at finite temperatures are the driving force for the transition. This is demonstrated in section III on the basis of Monte Carlo simulations of an Ising-type spin model, showing that a competition between AF Fe-Fe exchange-interactions and a nonmagnetic Rh state on the one hand and FM Fe-Rh interactions on the other is sufficient to explain the metamagnetic transition of FeRh. Inclusion of spatial degrees of freedom (section IV) by adding pair potentials and assuming a linear variation of the exchange parameter with the lattice parameter proves that this model explanation is in accordance with the ab initio data and leads furthermore to a nearly quantitative description of experimental details.

## II Ab initio total energy calculations

We have performed ab initio calculations by using the ASW-method Williams et al. (1979) and GGA. Relativistic effects are included in the scalar relativistic approximation. The basis wave functions of Fe and Rh atoms include spd and f states which are sufficient to obtain the correct magnetic behavior of $`\gamma `$-iron. We assume that the AF like spin structure can best be described by a spin-spiral (Ref. Uhl et al., 1994) with wave vector $`(0.5,0.5,0.5)`$ in units of $`2\pi /a`$. As first step, we optimized the volume of the AF-II state with equal ASA-radii for both types of atoms. For the resulting equilibrium volume, the total energy is then minimized with respect to the ratio of the ASA-spheres leading to $`r_{\mathrm{Rh}}/r_{\mathrm{Fe}}=1.15`$. In order to investigate the influence of the Rh moment in the FM phase we have also used the fixed spin moment method Uhl et al. (1994) to restrict $`m_{\mathrm{Rh}}`$ to zero.

Table 1 contains the calculated equilibrium properties for AF-II and FM phases which are in fair agreement with previous calculations. Compared to previous non-relativistic ASW LDA calculations (Moruzzi et al.Moruzzi and Marcus (1992a)), we obtain equilibrium volumes which are about 1-2$`\%`$ larger, being typical for GGA. In contrast to the cases of pure ironHerper et al. (1999) and $`\mathrm{FeNi}`$ Invar alloysBarbiellini et al. (1991), the magnetic energy differences are only slightly influenced by the gradient corrections. Also, no evidence was found that a noncollinear structure could be lower in energy than the previously found collinear AF and FM spin arrangements. With respect to the density of states at the Fermi level our results compare better with results of Szajek et al.Szajek and Morkowski (1994) than with results of Ref. Moruzzi and Marcus, 1992a. Since the latter results were obtained by using LDA ASW (with a spd basis), we conclude that the density of states at the Fermi level is very sensitive to computational details. The values obtained for the bulk moduli, $`B_{\mathrm{AF}}`$ and $`B_{\mathrm{FM}}`$, are lower compared to those of Ref. Moruzzi and Marcus, 1992a (where $`B_{\mathrm{AF}}>B_{\mathrm{FM}}`$ was obtained in contrast to $`B_{\mathrm{FM}}>B_{\mathrm{AF}}`$ in Ref. Szajek and Morkowski, 1994).

The calculated total energy curves are shown in Fig. 2. The large energetic difference between the usual FM phase and the hypothetic FM phase with zero rhodium moment implies that a finite Rh moment plays an important role for the stability of the FM phase. It seems unlikely that a magnetic field at the Rh sites induced by the surrounding iron atoms is responsible for the appearance of a Rh moment in the FM phase, as was proposed in previous discussions. However, a strong ferromagnetic exchange interaction between the Rh and the iron atoms that overrides an antiferromagnetic exchange between next nearest neighbor iron sites, would explain how the existence of a Rh moment can help to stabilize the FM phase. We have evidence as shown in the following sections, that a competition between a low lying nonmagnetic Rh state and another one with higher energy and finite moment (which can benefit from exchange with ferromagnetic iron neighbors) is the main reason for the metamagnetic transition.

## III A spin-analogous model

In order to confirm our hypothesis we have constructed a model being suitable for an examination of the metamagnetic transition at finite temperatures. The simplest way to do this is by means of Monte Carlo simulations with a localized spin model. To keep the model tractable, we neglect spin wave excitations and restrict ourselves to Ising spins. This is justified because the transition takes place between two ordered structures and both phases have collinear spin structures. For a description of the nonmagnetic Rh state in addition to the magnetic Fe and Rh states, we chose a spin-1 Ising model, where the spin variables can take the values $`S_i=1,0,+1`$. The spins are located on a bcc lattice with nearest and next nearest neighbor interactions. Depending on their positions, we distinguish between Fe or Rh sites, where each type occupies a simple cubic sublattice, corresponding to an ordered equiatomic alloy. The interaction parameters depend then on the type of sites involved. This situation can be described by the following Hamiltonian,

$$H=\underset{i}{\overset{}{}}D_iS_i^2\underset{<\mathrm{nn},\mathrm{nnn}>}{}J_{ik}S_iS_k.$$
(1)
Without the assumption of different types of atoms, Hamiltonian (1) is also known as the Blume-Capel model.Blume (1966); Capel (1966) The first term separates the nonmagnetic $`S_i=0`$ and the magnetic $`S_i=\pm 1`$ states. For Fe we choose a large positive value in order to suppress the $`S_i=0`$ state. For Rh we choose a negative value leading to a nonmagnetic ground state. The second term contains the exchange parameters $`J_{ik}`$ which depend only on the type of atoms located at the sites $`i`$ and $`k`$. In the case of ordered equiatomic FeRh, we have only three different parameters: $`J_{\mathrm{FeFe}}^{\mathrm{nnn}}`$, $`J_{\mathrm{FeRh}}^{\mathrm{nn}}`$ and $`J_{\mathrm{RhRh}}^{\mathrm{nnn}}`$. The first one is chosen to be negative in order to accomplish an AF ground state. The second is taken to be large and positive as outlined in the previous section. The third one, for the sake of simplicity, is set equal to zero. The choice for $`J_{\mathrm{FeFe}}^{\mathrm{nnn}}`$ is fixed by the Néel-temperature $`T_\mathrm{N}`$ of the AF phase assuming that no transition to a FM phase takes place. $`T_\mathrm{N}`$ can be determined from the P-T phase diagram by extrapolating the transition line between AF and PM phases (occurring at pressures larger than $`6\mathrm{GPa}`$) to zero pressure. $`J_{\mathrm{FeRh}}^{\mathrm{nn}}`$ and $`D_{\mathrm{Rh}}`$ have been chosen to yield realistic values for $`T_\mathrm{C}`$ and $`T_\mathrm{m}`$, respectively. The values of parameters used in the simulations are given in Table 2.

### III.1 Details of computation

The evaluation of thermodynamic properties of (1) is done on the basis of Monte Carlo simulations according to the Metropolis schemeMetropolis et al. (1953) using a sequential update. Interesting quantities like magnetization or magnetic moment are computed and summed up every $`10`$ to $`20`$ lattice sweeps, which ensures that the evaluated lattice configurations are sufficiently uncorrelated. Furthermore, we discard the first $`\mathrm{20\hspace{0.17em}000}`$ lattice sweeps in order to allow the system to reach thermal equilibrium before computing averages. In order to speed this up, we have also used the final configuration of the last run to initialize the simulation for the next temperature. Simulations which involve a phase transition are performed twice, with increasing and decreasing temperatures, in order to assure that thermal equilibrium has been reached.

The computed AF and FM phases are metastable, i. e. they are separated by a large energy barrier, which arises from the fact that in the transition states a considerable amount of FM domains have to be created in the AF phase and vice versa. So the standard algorithm is unlikely to overcome this barrier and as a result the metamagnetic transition might not be seen at all. Instead the phases have to be overheated or undercooled before transforming, which results in a large hysteresis or irreversible behavior which makes it difficult to obtain reliable information about the transition point. Therefore, it is necessary to modify the algorithm in order to allow for a direct jump to the other phase by circumventing the energy barrier with a global update step, where all spins will be updated at once. This algorithm ideally connects equilibrium configurations of the AF phase directly with equilibrium configurations of the FM phase while ensuring that the entropy difference between the states is correctly reproduced. This can be done by choosing a unique mapping between lattice configurations of the AF and FM phases, respectively. Or, more general, the selection probability of a specific target configuration must be the same as the selection probability of the previous start configuration in the backward direction.

Since it is a priori not clear how equilibrium configurations will look like at finite temperatures, we simply use an update scheme which connects the ground state configurations of both phases and thus works at least at low temperatures. In the vicinity of the transition temperature, this algorithm might not reproduce equilibrium states for the trial configuration, because nature and amount of excitations are presumably different in both phases. Since off-equilibrium means that most of the trial states are too high in energy, the probability that the trial state is by chance close enough to an equilibrium state to be accepted, decreases for larger system sizes.

In order to obtain a trial configuration, we divide the system into Rh and Fe sublattices. The Fe sublattice is again divided into two sublattices according to spin up and spin down positions (see Fig. 1). For each of these sublattices, it is decided randomly (with probability $`1/2`$), whether the corresponding sublattice is flipped as a whole. For the Rh sublattice another random number decides whether the spin values $`S_i\{1,0,+1\}`$ are rotated clockwise or counter-clockwise (i. e. $`1`$ becomes $`0`$, $`0`$ becomes $`+1`$, $`+1`$ becomes $`1`$, or vice versa). So, each trial configuration is chosen with the same probability. Afterwards the energy difference $`\mathrm{\Delta }H`$ between the present and the previous configuration is computed and the new configuration is accepted with probability $`\mathrm{min}(\mathrm{exp}(\mathrm{\Delta }H/k_\mathrm{B}T),1)`$ assuring detailed balance. This global update step is performed after each complete lattice sweep.

We performed for each temperature between $`\mathrm{100\hspace{0.17em}000}`$ and $`\mathrm{1\hspace{0.17em}000\hspace{0.17em}000}`$ (around $`T_\mathrm{m}`$) lattice sweeps for different linear system sizes $`L=6`$ to $`16`$. However, our global update scheme does only show a metamagnetic transition within reasonable simulation times up to a system size of $`L=10`$. Looking for a different way to determine the transition point which would also work for larger system sizes, we estimate the free energy by integrating the specific heat,

$$F(T)=E(T)T\left(S(0)+_0^T\frac{C}{T}𝑑T\right).$$
(2)
In the computation we fit the simulated results for the specific heat divided by temperature, $`C/T`$, and the internal energy $`E`$ with 10th order polynomials, which can be integrated analytically. For the system characterized by the Hamiltonian (1), we neglect the entropy contribution at zero temperature $`S(0)`$, since the ground state spin structure in both phases is nondegenerate except for systems with spin inversion symmetry. Furthermore, concerning the estimation of the free energy finite size effects are not expected to affect the results, since phase transitions are not encountered during these simulations. For the rest of the calculations a comparison of the results for smaller and larger systems sizes (as can be seen in the upcoming figures) reveals that the main issues of this paper are also not affected by the restricted system sizes.

### III.2 Computational results

The order parameter of the AF phase is the staggered magnetization for the AF-II spin structure, i. e. the sum (of the absolute values) of the staggered magnetizations of the two simple cubic sublattices that constitute the bcc lattice structure. The staggered magnetization of a simple cubic (NaCl) lattice is defined as sum of spins multiplied with a sign which alternates depending on whether the corresponding spin occupies a Na or a Cl position. The order parameter of the FM phase is given by the magnetization of the lattice. The variation of both order parameters with temperature is shown in Fig. 3 for a system of size $`L=10`$. At low temperatures the staggered magnetization approaches a maximum value of $`0.5`$ due to the fact that the Rh sublattice has no moment (Fig. 4) and therefore does not contribute to the sum. At $`T_\mathrm{m}=268`$K the staggered magnetization abruptly drops to zero, whereas magnetic moment and magnetization increase close to their saturation values. Above $`T_\mathrm{C}=610`$K ferromagnetism breaks down and the system becomes paramagnetic. At the same time, the average moment of the Rh atoms falls down to a value of $`0.2\mu _\mathrm{B}`$. This is also the reason, why the decrease in the magnetization appears unusually sharp. From the data, a phase transition of second as well as of first order seems to be possible. Accordingly, it is known that for negative values of $`D`$ the Blume-Capel model can show both kinds of phase transitions separated by a tricritical point in the $`T`$-$`D`$-plane.Blume et al. (1971); Saul et al. (1974); Jain and Landau (1980) Hence in order to safely determine the nature of this phase transition further calculations are necessary.

A close look at the specific heat (Fig. 5) helps to explain the occurrence of a metamagnetic transition. At first sight it seems as if the specific heats of FM and AF phases do not differ very much below $`350`$K. Above this temperature the specific heat of the AF phase increases more rapidly with temperature until around $`570`$K the overheated AF phase is not stable anymore in the simulations and transforms to the FM phase. However, the inset shows that starting around $`100`$K the specific heat of the FM phase is enhanced compared to the AF heat. We explain this enhancement by a Schottky-type anomaly which adds up to the excitations from spin flips. Schottky anomalies are observed in systems with two levels separated by a small energy barrier. In fact a crossover from magnetic to nonmagnetic Rh atoms is conform with this picture. In the FM phase the Rh atoms have a moment being ferromagnetically aligned to the Fe moments, because the loss of the moment would correspond to an energy loss of eight times the exchange constant. This amount is diminished by the energy gain due to the $`D_{\mathrm{Rh}}`$ term, which is smaller than $`8\times J_{\mathrm{FeRh}}^{\mathrm{nn}}`$ but larger than $`4\times J_{\mathrm{FeRh}}^{\mathrm{nn}}`$, since the ground state is otherwise not antiferromagnetic. Increasing fluctuations of the magnetization of the Fe sublattice cause this energy difference to decrease and finally lead to a breakdown of the average Rh moment at $`T_\mathrm{C}`$. In the AF phase magnetic Rh atoms can be excited at the expense of the energy $`D_{\mathrm{Rh}}`$. There is, however, no gain in energy due to the exchange interaction, since the contributions from the AF Fe atoms cancel at the Rh site. This corresponds to a much larger energy difference between magnetic and nonmagnetic Rh states compared to the FM phase and does not lead to an appreciable contribution to the specific heat.

This view is further supported by a comparison of the free energy of both phases shown in Fig. 6: The curves intersect at $`T266`$K, which is in excellent agreement with the simulated $`T_\mathrm{m}=268`$K. This accuracy could be achieved, because the data for the specific heat obtained from the simulations have only little spread and the fitted curves interpolate the data points perfectly. From the difference between the internal energies of both phases at the transition point, the entropy jump at $`T_\mathrm{m}`$ is determined to $`\mathrm{\Delta }S(T_\mathrm{m})=5.2\mathrm{J}\mathrm{Kg}^1\mathrm{K}^1`$, which is only $`30`$ to $`45\%`$ of the experimental values (but of the right order of magnitude). But, one has to bear in mind that with the choice of an Ising model we neglected the possibility of non-collinear moments. For the Rh moments this can only occur in the FM phase, and would therefore contribute to $`\mathrm{\Delta }S(T_\mathrm{m})`$. A second point is that the weight of the $`S_i=0`$ state is the same as for each magnetic state $`S_i=\pm 1`$. This is a natural choice for the spin-1 Ising model; but for the real system this is somehow arbitrary, because we have no information about the electronic origin of both Rh states.

## IV An extended model

Since the simple spin Hamiltonian (1) can reproduce a metamagnetic transition as observed in $`\alpha `$-FeRh, it remains an interesting question whether other outstanding properties of this alloy, as the large volume increase at $`T_\mathrm{m}`$, can also be explained. Furthermore, it has still to be proven whether our spin analogy is a good approximation to the ab initio results in the sense that it has comparable low temperature properties. In order to check this we have to extend (1) for a description of elastic and magneto-volume properties:

$`H`$ $`=`$ $`{\displaystyle \underset{i}{\overset{}{}}}D_iS_i^2{\displaystyle \underset{<\mathrm{nn},\mathrm{nnn}>}{}}J_{ik}(r_{ik})S_iS_k`$ (3)
$`+{\displaystyle \underset{<\mathrm{nn}>}{}}V_{\mathrm{nn}}(r_{ik})+{\displaystyle \underset{<\mathrm{nnn}>}{}}V_{\mathrm{nnn}}(r_{ik}).`$
For $`V(r_{ik})`$ we use simple pair potentials of the Lennard-Jones type:

$$V(r_{ik})=4ϵ\left[\left(\frac{\sigma }{r_{ik}}\right)^{12}\left(\frac{\sigma }{r_{ik}}\right)^6\right].$$
(4)
Since in general the lattice structure of a Lennard-Jones system is closely packed, two different pair potentials for nearest and next-nearest neighbors have to be used in order to stabilize the bcc structure. The potentials, however, do neither distinguish between different atom types nor between the different spin states, as has been done in previous simulations of related materials like Fe-Ni Invar or Y(MnxAl1-x)2.Gruner et al. (1998); Gruner and Entel (2002) The use of Lennard-Jones potentials is far from being optimum for metals, but has numerical advantages that enable us to speed up the calculations substantially. It is then sufficient to choose the parameters so that basic elastic properties like low temperature lattice constant, bulk modulus or thermal expansion are reproduced.

Another change compared to (1) is that the exchange parameter is now taken to be a function of interatomic distance. For simplicity, we assume a linear distance dependence:

$$J_{ik}(r_{ik})=J_{ik}+\frac{J_{ik}}{r}r_{ik}.$$
(5)
The values for $`J_{ik}`$ which we take in the extended calculations are roughly the same as in Table 2. The relation between the derivatives of $`J_{\mathrm{FeRh}}^{\mathrm{nn}}`$ and $`J_{\mathrm{FeFe}}^{\mathrm{nnn}}`$ with respect to the interatomic distance was determined by the relation of the pressure derivatives of the Curie and Néel temperatures, respectively, which have been obtained from the experimental phase diagram by assuming a linear dependence between $`J_{ik}/r`$ and $`T_{\mathrm{C},\mathrm{N}}/p`$. The absolute values have been adapted to reproduce the volume jump at $`T_\mathrm{m}`$. An exchange interaction of this form has also been used to describe magneto-volume effects in Fe-Ni Invar.Grossmann and Rancourt (1996); Rancourt and Dang (1996); Lagarec and Rancourt (2000) In this case, the derivative of the exchange constant $`J_{ik}/r`$ was a factor of $`2`$ to $`20`$ larger than in the present work. The values of parameters are summarized in Table 3.

### IV.1 Details of computation

For the evaluation of Hamiltonian (3) we use a textbook isothermal-isobaric Monte Carlo method (e. g. Ref. Allen and Tildesley, 1987) consisting of alternating spin and position updates for each atom and a global volume update step after finishing each lattice sweep. This algorithm has been used by the the authors in previous calculations and is explained in detail in the corresponding references.Gruner et al. (1998); Gruner and Entel (1998, 2002) In order to simulate the metamagnetic transition, an additional global spin update step has to be introduced. We use the same algorithm as described in the last section for the spin system in connection with a simultaneous volume adaption. Since the latter reduces the acceptance probability considerably, the global update scheme has to be repeated several thousands of times. It is not practicable to use a new spin configuration for each trial step, because the evaluation of the energy is comparatively time consuming. On the other hand, the new energy after solely rescaling the volume can be calculated very quickly, since due to the use of Lennard-Jones potentials the energy can be written as a function of integral powers of the lattice parameter. Therefore, we choose a trial spin configuration as described before and compute the new energy. Then we attempt a previously fixed number $`N`$ of Metropolis steps, each with a newly chosen volume. If one step is accepted, we continue with the original spin configuration as trial system (and so on) until $`N`$ steps have been made. Since the number of trial steps $`N`$ has been previously fixed, detailed balance is still valid.

We have performed simulations with system sizes of $`L=6`$ to $`12`$. A direct metamagnetic transition could only be seen for system sizes up to $`L=8`$. The simulation time ranged from $`\mathrm{120\hspace{0.17em}000}`$ up to $`\mathrm{1\hspace{0.17em}000\hspace{0.17em}000}`$ lattice sweeps around $`T_\mathrm{m}`$ with values from $`1000`$ to $`\mathrm{10\hspace{0.17em}000}`$ for $`N`$. As before we estimated Gibbs’ free energy for zero pressure by integrating the specific heat according to Eq. (2).

### IV.2 Computational results

For a comparison of our model properties with the results of ab initio calculations, we calculated in isochoric simulations the energy as a function of the volume at low temperatures for different fixed spin structures (Fig. 7). Here, we cooled a system of size $`L=16`$ exponentially down from $`100`$K down to $`4`$K. We find good agreement with the results of Fig. 2 in the sense that the order of magnetic phases is similar. Although the ab initio total energy differences are somewhat larger, this indicates that our model Hamiltonian (3) is a qualitatively correct description of the mechanisms leading to a metamagnetic transition in FeRh.

As in the simple model without volume-dependent terms, we find an abrupt increase of the magnetization in combination with a discontinuous decrease for the staggered magnetization with increasing temperature (Fig. 8). Consequently the mean moment at the Rh sites (Fig. 9) also raises sharply around the metamagnetic transition temperature of $`T_\mathrm{m}=322`$K. Around the Curie temperature $`T=720`$K the magnetization decreases more smoothly than for the simple model. This may be due to the enhancement of the effective exchange parameter given by Eq. (5) which is caused by the lattice expansion.

As expected from the low temperature calculations, the volume of the AF phase is smaller than the volume of the FM phase throughout the stability range (Fig. 10). For the freely fluctuating system a volume jump of $`0.8\%`$ occurs at $`T_\mathrm{m}`$. Below $`T_\mathrm{C}`$ the volume expansion is reduced in the FM phase, which is in qualitative agreement with experimentZakharov et al. (1964); Aigarabel et al. (1996).

The bulk modulus $`B_\mathrm{S}`$ of the AF phase is about $`6\%`$ larger than the FM bulk modulus throughout the stability range (Fig. 11). Around $`T_\mathrm{m}`$ and $`T_\mathrm{C}`$ we find a considerable weakening of the material which is caused by the magneto-volume anomalies. The values for $`B_\mathrm{S}`$ at low temperatures can be estimated more accurately by using a polynomial fit to the $`E`$-$`V`$ curves in Fig. 7. We find $`B_\mathrm{S}=231`$GPa for the AF phase and $`B_\mathrm{S}=217`$GPa for the FM phase. These values are in good agreement with results of isobaric calculations shown in Fig. 11. Compared with the results of ab initio calculations (Table 1) the absolute values are too large, whereas these calculations do not give a unanimous prediction for the sign and the magnitude of the difference $`\mathrm{\Delta }B_\mathrm{S}`$ (with respect to the magnetic structures). Experimental measurementsZakharov et al. (1964) of Young modulus suggest that the bulk modulus of the AF phase should in fact be lower than in the FM phase. The modulus of the AF phase, however, has only been estimated in the vicinity of $`T_\mathrm{m}`$, where a weakening of the material (as a precursor of the transition) is present (as in our simulations).

The specific heats obtained for the extended model (Fig. 12) resemble the findings for the spin-only Hamiltonian (1), except that now additional contributions due to atomic displacements are included. This shows that the proposed explanation of the metamagnetic transition is still valid. Since the Hamiltonian (3) does not contain any kinetic terms, the low temperature value of the specific heat is only half of the Dulong-Petit limit of $`3R`$, where $`R`$ is the kinetic gas constant.

### IV.3 Estimation of Gibbs free energy

The calculation of Gibbs free energy $`G`$ for the model with classical motion (3) is more complicated than in the spin-only case (1). First of all absolute values of $`G`$ cannot be given, since the specific heat at zero temperature is finite and the resulting entropy would diverge. But since the limit for $`T`$$``$$`0`$ of $`C_\mathrm{p}`$ is $`1.5R`$ and hence the same for all magnetic structures, differences of free energies can be computed. In contrast to the previous section, the entropy contribution at zero temperature $`\mathrm{\Delta }S(0)`$ must be considered. As before a contribution from the magnetic system can be neglected, while we have to account for the differences in the elastic properties of FM and the AF phases. Since we only need the entropy at zero temperature, we choose an ensemble of harmonic oscillators as an approximation for our spatial degrees of freedom, neglecting the anharmonicity of the potentials and the coupling of the oscillators. Differentiating the free energy obtained from the logarithm of the partition function leads then to a simple expression for the entropy difference,

$$\mathrm{\Delta }S=\frac{3}{2}Nk_\mathrm{B}\mathrm{ln}\left(k_{\mathrm{AF}}/k_{\mathrm{FM}}\right),$$
(6)
where $`k_{\mathrm{AF}}`$ and $`k_{\mathrm{FM}}`$ are the force constants of the harmonic potentials which are estimated from the curvature of the ground state energy (versus lattice parameter) curve. For the entropy difference we obtain then: $`\mathrm{\Delta }S(0)=9.17\mathrm{J}\mathrm{kg}^1\mathrm{K}^1`$. Taking this into account, we achieve again a rather good value for $`T_\mathrm{m}`$. From the differences of internal energies at $`T_\mathrm{m}`$ we obtain for the entropy jump at the transition $`\mathrm{\Delta }S(T_\mathrm{m})=15.9\mathrm{J}\mathrm{kg}^1\mathrm{K}^1`$, which is within the range of experimental results. Comparison of the calculated free energy with experimental values for $`\mathrm{\Delta }G`$ shows excellent agreement. The experimental values have been obtained by a graphical integration of the magnetic field expressed as a function of the measured magnetization.Ponomarev (1972) Extrapolation of experimental values to zero temperature shows that the energy difference $`\mathrm{\Delta }E(0)`$ between the AF and FM phases is much better described by our model parameters than by ab initio results, since the latter show that $`\mathrm{\Delta }E`$ is one order of magnitude too large. This discrepancy has already been noticed by Moruzzi and MarcusMoruzzi and Marcus (1992b) who relate this to the omission of zero point energy corrections in their total energy calculations.

## V Conclusions

We propose on the basis of new ab initio results a mechanism for the metamagnetic transition in FeRh at finite temperatures. In contrast to previous explanations, our model does not rely on a large difference between the low temperature specific heat constants of both phases. These are expected to be sensitive to external influences as both experimental measurements and band structure calculations suggest, so that it seems implausible that a constant contribution of the given magnitude might survive up to room temperature. Instead, we propose that the existence of two magnetic states of Rh atoms connected with competing FM Fe-Rh and AF Fe-Fe exchange interactions are at the origin of the metamagnetic transition. The magneto-volume effects can simply be explained on the basis of distance dependent exchange parameters. The applicability of this mechanism to the Fe-Rh problem has been verified by Monte Carlo model calculations, showing that a metamagnetic transition of the desired kind does in fact occur and, by extending the model, magneto-volume effects and other experimental properties can be sufficiently well described.

As we have pointed out, our explanation is in agreement with existing experimental data. However, a further check of our model would be a comparison with (non-existing) specific heat data from above $`T_\mathrm{C}`$ down to very low temperatures for both, the AF and FM phases. From this one could then estimate the magnetic contribution by subtracting the lattice part within the Debye approximation, the electronic part and the contribution by the anharmonicity of the potentials.Bendick et al. (1978); Rellinghaus et al. (1995) So far, the specific heat has only been determined in the range from $`100`$K to $`500`$K and only for the nearly stoichiometric Rh-rich alloy.Richardson et al. (1973); Annaorazov et al. (1992) Systematic measurements on the Fe-rich side with a FM ground state and measurements under pressure, suppressing the metamagnetic transition would yield information whether a Schottky-type excitation plays an important role, which should show up around $`200`$K in the magnetic contribution to the specific heat of ferromagnetic samples.

Monte Carlo simulations with applied pressure are left for future work, since with increasing pressure and hence increasing $`T_\mathrm{m}`$ a reliable estimation of the transition temperature is rather difficult and requires an improvement of the global MC step. First tests, however, showed that the location of the phase boundaries under pressure is in sufficient agreement with experimental data for pressures below $`20`$kbar and above $`40`$kbar, where the metamagnetic transition is completely suppressed. The tricritical point, if one exists, should be located somewhere between $`20`$kbar and $`40`$kbar for the parameters used here.

###### Acknowledgements.

This work has been supported by the DFG (Deutsche Forschungsgemeinschaft) through the SFB (Sonderforschungsbereich) 445 and the Graduate College Structure and Dynamics of Heterogeneous Systems.


text:
# Isotropization in Brane Gas Cosmology

## I Introduction

The Standard Big Bang (SBB) cosmology has become an extremely successful model that has been well tested by experiment. However, the model is incomplete. The underlying theory of classical general relativity and the description of matter as an ideal gas breaks down at the high temperatures of the early Universe, and the solutions of the theory in fact have an initial singularity. Moreover, SBB does not address many important cosmological questions such as the observed homogeneity, spatial flatness, and the origin of structure in the Universe. Cosmological inflation (see e.g.  for textbook treatises and  for shorter reviews) builds on SBB cosmology providing a solution to some of these issues, but it (at least in the context of scalar field-driven inflation) suffers from the same initial singularity problem  and other conceptual problems , which indicate that inflation cannot be the complete story of early Universe cosmology.

In recent years, many models motivated by string theory and M-theory have emerged as possible solutions to the outstanding problems of early Universe and inflationary cosmology (see e.g.  for recent but incomplete reviews). Beginning with the work on Pre-Big-Bang cosmology  it was realized that a dynamical dilaton should play an important role in the very early Universe. More recently, models have become prominent in which our Universe consists of a 3-brane embedded in a higher dimensional bulk space, with the standard model constrained to live on the brane . Although these models can resolve a number of issues, such as the hierarchy problem, they introduce several other difficulties in the process. For example, large extra dimensions should be explained by classical general relativity, and it has been shown this results in problems stabilizing the brane . More importantly, in most of these models the six “extra” spatial dimensions are taken to be compactified, a priori, with no explanation provided for how this could come about dynamically. Although it appears to be an important concern for the naturalness of these models, this issue is rarely discussed in the literature.

An alternative approach to string/M-theory cosmology is the string gas or BV scenario. This model began with works  in which the effects of string gases on the cosmological evolution of the low energy effective string theory background geometry including the dilaton were explored. The most important result to emerge from these works is a dynamical mechanism, tied to the existence of string winding modes, which yields a nonsingular cosmology and may explain why at most three spatial dimensions can become large if the initial state is chosen to correspond to a Universe which is small in all spatial directions.

This work has been generalized to include the cosmological effects of p-brane gases and leads to the current model of Brane Gas Cosmology (BGC) . In BGC, the Universe starts analogous to the SBB picture, i.e. hot, dense, and with all fundamental degrees of freedom in thermal equilibrium. The Universe is assumed to be toroidal in all nine spatial dimensions and filled with a p-brane gas. The assumption of toroidal geometry of the background space leads to the existence of string winding modes, since the background space admits cycles on which branes of the relevant dimensionalities (in particular one branes) can wrap. This wrapping is associated with a winding energy which - in the context of dilaton gravity - acts as a confining potential for the scale factor preventing further expansion of the spatial dimensions. Also associated with the brane are oscillatory modes described by scalar fields and momentum modes which correspond to the center of mass motion of the brane. The momentum modes are related by T-duality to the winding modes and this duality results in the non-singular behavior of the model. In order for dimensions to decompactify, p-brane winding modes must annihilate with anti-winding modes and it is argued that this only occurs in a maximum of $`2p+1`$ dimensions . Since strings $`(p=1)`$ are the lowest dimensional objects that admit winding modes, since they are the lightest of all winding modes and hence fall out of equilibrium later than other winding modes, it follows that the number of large space-time dimensions can be at most (3+1).

In the context of the background equations of dilaton gravity, the winding modes yield a confining potential for the scale factor which also gives rise to a period of cosmological loitering (expansion rate near zero) for the three large spatial dimensions. This is due to the time needed for winding modes to annihilate and produce closed strings or loops . This is of great interest since loitering can explain the horizon and relic problems of standard cosmology without resorting to inflation However, to obtain a solution of the flatness and entropy problems, a phase of inflation following the decoupling of the three large spatial dimensions may be required. It was also shown in  that by considering loop production at late times BGC naturally evolves into the SBB, with a $`3+1`$ dimensional, radiation dominated Universe.

There are important issues that remain to be addressed within BGC. The fact that toroidal geometry was assumed for the background space was used for the existence and topological stability of winding modes. However, K3 or Calabi-Yau manifolds are more realistic choices for backgrounds within string theory and they do not admit 1-cycles (necessary to have topologically stable winding modes). Promising results have recently appeared which indicate that the conclusions of BGC extend to a much wider class of spatial background, including backgrounds which are K3 or Calabi-Yau manifolds .

Perhaps the main issue to be addressed in BGC is that of spatial inhomogeneities. That is, we would naturally expect fluxes and p-brane sources to lead to the possibility of catastrophic instabilities of spatial fluctuation modes. Although we do not address this issue here, we plan to study the role of inhomogeneities in followup work.

Other important issues for BGC are the questions of stabilization of the six small extra dimensions and isotropization of the three dimensions that grow large. Although these topics may seem to be unrelated, they both can be addressed by generalizing the BGC scenario to the anisotropic case. This paper will concentrate on the isotropization of the three large dimensions, but the generalization of BGC to the anisotropic case achieved in this paper will be valuable to address the issue of stabilization in later work.

## II Brane Gas Cosmology

We begin with a brief review of BGC, for more details the reader is referred to . Consider compactification of 11-dimensional M-theory on $`S^1`$, which yields type II-A string theory in 10-dimensions. The fundamental degrees of freedom in the theory are 0-branes, strings, 2-branes, 5-branes, 6-branes, and 8-branes. The low energy effective action of the theory is that of supersymmetrized dilaton gravity,

$$S_{bulk}=\frac{1}{2\kappa ^2}d^{10}x\sqrt{G}e^{2\varphi }\left[R+4G^{\mu \nu }_\mu \varphi _\nu \varphi \frac{1}{12}H_{\mu \nu \alpha }H^{\mu \nu \alpha }\right],$$
(1)
where $`G`$ is the determinant of the background metric $`G_{\mu \nu }`$, $`\varphi `$ is the dilaton given by the radius of the $`S^1`$ compactification, $`H`$ denotes the field strength corresponding to the bulk antisymmetric tensor field $`B_{\mu \nu }`$, and $`\kappa `$ is determined by the ten-dimensional Newton constant.

Fluctuations of each of the $`p`$-branes are described by the Dirac-Born-Infeld (DBI) action  and are coupled to the ten-dimensional action via delta function sources. The DBI action is

$$S_p=T_pd^{p+1}\zeta e^\varphi \sqrt{det(g_{mn}+b_{mn}+2\pi \alpha ^{}F_{mn})},$$
(2)
where $`T_p`$ is the tension of the brane, $`g_{mn}`$ is the induced metric on the brane, $`b_{mn}`$ is the induced antisymmetric tensor field, and $`F_{mn}`$ is the field strength tensor of gauge fields $`A_m`$ living on the brane. The constant $`\alpha ^{}l_{st}^2`$ is given by the string length scale $`l_{st}`$.

In this paper we will concentrate on the role of fundamental strings ($`p=1`$), since strings play the critical role in the dynamics . We will also ignore the effects of fluxes and string oscillatory modes since we are working in the low energy regime. Our basic approach is to consider the string gas as a matter source for the dilaton-gravity background. The string winding modes play a key role. The main point is that we want to find the first departures from the standard FRW picture resulting from considering stringy effects. This should result from the winding modes of the strings and the presence of the dilaton in the equations of motion.

In this paper, we generalize the equations of BGC to the anisotropic case. Thus, we take the metric to be of the form

$$ds^2=dt^2\underset{i=1}{\overset{D}{}}e^{2\lambda _i}(dx_i)^2,$$
(3)
where the label $`i`$ runs over the $`D`$ spatial indices, $`x_i`$ are the comoving spatial coordinates, $`t`$ is physical time, and the scale factor in the i’th direction is $`log(\lambda _i)`$. We again stress that inhomogenities are of vital importance, however we will leave their investigation to future work and take the string gas to be homogeneous.

By varying the total action we obtain the following equations,

$`{\displaystyle \underset{i=1}{\overset{D}{}}}\dot{\lambda }_i+\dot{\phi }^2=e^\phi E`$ (4)
$`{\displaystyle \underset{i=1}{\overset{D}{}}}\dot{\lambda }_i+\ddot{\phi }={\displaystyle \frac{1}{2}}e^\phi E`$ (5)
$`\ddot{\lambda }_i\dot{\phi }\dot{\lambda }_i={\displaystyle \frac{1}{2}}e^\phi P_i,`$ (6)
where $`E`$ is the total energy of the strings and $`P_i`$ is the pressure in the $`ith`$ direction, and we have introduced a “shifted” dilaton field

$$\phi =\mathrm{\hspace{0.17em}2}\varphi \underset{i=1}{\overset{D}{}}\lambda _i.$$
(7)
We can see immediately from the first equation that $`\phi `$ can not change sign. To insure that the low energy approximation remains valid we choose $`\phi <0`$ and $`\dot{\phi }>0`$. This causes the dilaton to have a damping effect in the equations of motion Alternatively we could consider $`\phi >0`$ and one might anticipate that this could lead to inflation. However this is a false conclusion, since the presence of winding modes act as a confining potential which can be seen from the negative pressure term in equation (6).. Another important observation is that these equations reduce to the standard FRW equations for fixed dilaton as was discussed in . This is comforting because at late times it is expected that the dilaton will gain a mass, perhaps associated with supersymmetry breaking and as a result we naturally regain the usual SBB cosmology. The string winding states will appear at late times as solitons .

### A Winding Modes

The energy and pressure terms follow from varying the brane action (2). The energy associated with the winding modes in the $`ith`$ direction is given by (taking the winding strings to be straight for simplicity)

$$E_i^w=\mu N_i(t)a_i(t)=\mu N_i(t)e^{\lambda _i},$$
(8)
where $`N_i`$ is the number of winding modes in the ith direction and $`\mu `$ is the mass per unit length of the string multiplied by the initial spatial dimension. It follows that the total energy is given by,

$$E_T^w=\underset{i=1}{\overset{D}{}}\mu N_i(t)e^{\lambda _i}.$$
(9)
In the same approximation that the winding strings are straight, the corresponding pressure terms are given by ,

$`P_i^w`$ $`=`$ $`\mu N_i(t)e^{\lambda _i}`$ (10)
$`P_j^w`$ $`=`$ $`0(ji).`$ (11)
It follows then by inserting (10) into (6) that winding strings give rise to a confining potential for the scale factor, and hence prevent expansion. In order for a spatial dimension to be able to become large, the winding modes in that direction must be able to annihilate.

### B Annihilation and Loop Production

As argued earlier, the fact that strings are likely to annihilate in a maximum of three space dimensions leads to three of the initial nine spatial dimensions growing large. As the background continues to expand the three dimensional space will be filled with loops resulting from the annihilation of the winding states. The corresponding energy loss in winding states and the energy and pressure of the creation of loops must be taken into account when considering the dynamics. The strings in the expanding space will become of macroscopic size and the required pressure and energy terms will be analogous to that of a cosmic string network .

We know that the energy in winding modes is given by (9). Considering the time rate of change of this energy we find,

$$\dot{E}_T^w=\underset{i=1}{\overset{D}{}}\mu \dot{N}_i(t)e^{\lambda _i}+\underset{i=1}{\overset{D}{}}\mu N_i(t)e^{\lambda _i}\dot{\lambda }_i.$$
(12)
The first term on the right hand side of this equation corresponds to energy loss into string loops or radiation, the second term corresponds to the gain in total energy of the strings due to the stretching by the expansion of space.

The energy in loops can be written as

$$E^{loops}=g(t)V_0,$$
(13)
where $`g(t)`$ is the energy density in string loops (plus radiation) per initial comoving volume (and is a constant if no loop production or energy loss from winding strings into radiation occurs), and $`V_0=\mathrm{exp}(\lambda _{0i})`$ is the initial volume, $`\lambda _{0i}`$ being the initial values of the logarithms of the scale factors. Since sting loops are likely to be produced relativistically, we will use the equation of state of radiation to describe them. This approximation also allows us to treat the string loops and other radiation together. This approximation could easily be relaxed without changing our basic conclusions. Thus, we use the equation of state of relativistic radiation $`p=(1/3)\rho `$ ($`\rho `$ and $`p`$ denoting energy density and pressure, respectively) which implies,

$$P^{loops}=\frac{1}{3}g(t)V_0.$$
(14)
To find how $`g(t)`$ evolves we equate the loss in winding energy due to energy transfer (the first term on the right hand side of 12) with the change in loop energy,

$$\underset{i=1}{\overset{D}{}}\mu \dot{N}_i(t)e^{\lambda _i}=\dot{g}(t)V_0.$$
(15)
Thus, the loop production is determined by the change in the number of winding modes, as expected:

$$\dot{g}(t)=V_0^1\underset{i=1}{\overset{D}{}}\mu \dot{N}_i(t)e^{\lambda _i}.$$
(16)
It remains to determine the rate of winding mode annihilation, but before doing so let us consider the specific case of a 2+1 anisotropic Universe.

## III Anisotropic Generalization

As three of the nine spatial dimensions grow large there is no a priori reason to expect that this should happen in an isotropic manner. Moreover, since the dimensions expanding correspond to the annihilation of winding modes due to string intersections, we might expect this process to occur at differing rates leading to anisotropic dimensions. However, we might also expect that as the winding modes annihilate in one dimension and that dimension begins to expand faster there are then fewer winding modes left to annihilate. In this way the remaining dimensions are given an opportunity to isotropize. To explore if this is indeed the case, let us consider the case of a 2+1 anisotropic background, whose metric we write in the form

$$ds^2=dt^2e^{2\lambda }(dx^2+dy^2)e^{2\nu }dz^2.$$
(17)
The scale factor corresponding to $`\lambda `$ is denoted $`a(t)`$, the one corresponding to $`\nu `$ is denoted $`b(t)`$.

The equations of motion (4) then become,

$`2\dot{\lambda }^2\dot{\nu }^2+\dot{\phi }^2=e^\phi E`$ (18)
$`2\dot{\lambda }^2\dot{\nu }^2+\ddot{\phi }={\displaystyle \frac{1}{2}}e^\phi E`$ (19)
$`\ddot{\lambda }\dot{\phi }\dot{\lambda }={\displaystyle \frac{1}{2}}e^\phi P_\lambda `$ (20)
$`\ddot{\nu }\dot{\phi }\dot{\nu }={\displaystyle \frac{1}{2}}e^\phi P_\nu .`$ (21)
The energy and pressure terms are given by

$`E=E_T^w+E^{loops}`$ $`=`$ $`2\mu N(t)e^\lambda +\mu M(t)e^\nu +g(t)e^{2\lambda _0+\nu _0},`$ (22)
$`P_\lambda `$ $`=`$ $`\mu N(t)e^\lambda +{\displaystyle \frac{1}{3}}g(t)e^{2\lambda _0+\nu _0},`$ (23)
$`P_\nu `$ $`=`$ $`\mu M(t)e^\nu +{\displaystyle \frac{1}{3}}g(t)e^{2\lambda _0+\nu _0},`$ (24)
where $`N(t)`$ and $`M(t)`$ are the number of winding modes in the $`\lambda `$ and $`\nu `$ directions, respectively.

We now consider the effect of loop production as a result of winding mode annihilation. To simplify the analysis, we will assume that winding modes in $`\lambda `$ and $`\nu `$ directions form two separate, noninteracting gases. We expect that this approximation will reduce the isotropization of winding modes, and hence that including the interactions we omit will lead to accelerated isotropization. This issue is being studied currently .

Winding mode annihilation results from string intersections, which depends  on the square of the number of modes, inversely on the cross-sectional area available for the interaction and is proportional to the Hubble length of the string. Since causality plays a role in the interaction, we must scale our dimensions by the Hubble time $`t`$. Hence,

$$\dot{N}(t)=\stackrel{~}{c}_NN^2(t)\left(\frac{t^2}{Area}\right)t^1=\stackrel{~}{c}_NN^2(t)\left(\frac{t^2}{a(t)b(t)}\right)t^1=c_NtN^2(t)e^{\lambda \nu },$$
(25)
where $`\stackrel{~}{c}_N`$ is a dimensionless constant, and $`c_N`$ is this same constant rescaled by the basic length dimensions related to $`a`$ and $`b`$. Similarly, for $`M(t)`$ we have

$$\dot{M}(t)=c_MtM^2(t)e^{2\lambda }.$$
(26)
Using these expressions in (16) gives us the evolution of loop production,

$$\dot{g}(t)=\mu te^{2\lambda _0\nu }\left(2c_NN^2(t)e^\nu +c_MM^2(t)e^{2\lambda +\nu }\right).$$
(27)
For computational simplicity we define $`l(t)\dot{\lambda }`$, $`q(t)\dot{\nu }`$, and $`f(t)\dot{\phi }`$, which leaves us with the following set of first order ordinary differential equations,

$`2l^2q^2+f^2=e^\phi \left(2\mu N(t)e^\lambda +\mu M(t)e^\nu +g(t)e^{2\lambda _0+\nu _0}\right)`$ (28)
$`\dot{f}=2l^2+q^2+{\displaystyle \frac{1}{2}}e^\phi \left(2\mu N(t)e^\lambda +\mu M(t)e^\nu +g(t)e^{2\lambda _0+\nu _0}\right)`$ (29)
$`\dot{l}=l\phi +{\displaystyle \frac{1}{2}}e^\phi \left(\mu N(t)e^\lambda +{\displaystyle \frac{1}{3}}g(t)e^{2\lambda _0+\nu _0}\right)`$ (30)
$`\dot{q}=q\phi +{\displaystyle \frac{1}{2}}e^\phi \left(\mu M(t)e^\nu +{\displaystyle \frac{1}{3}}g(t)e^{2\lambda _0+\nu _0}\right)`$ (31)
$`\dot{N}(t)=c_NtN^2(t)e^{\lambda \nu }`$ (32)
$`\dot{M}(t)=c_MtM^2(t)e^{2\lambda }`$ (33)
$`\dot{g}(t)=\mu te^{2\lambda _0\nu }\left(2c_NN^2(t)e^\nu +c_MM^2(t)e^{2\lambda +\nu }\right)V_0^1`$ (34)
$`l(t)=\dot{\lambda }`$ (35)
$`q(t)=\dot{\nu }`$ (36)
$`f(t)=\dot{\phi }.`$ (37)
The system (28-37) is over-determined (there are more equations than unknowns). Therefore we take (28) as a constraint on the initial conditions. We find from numerical analysis that the winding modes $`N(t)`$ and $`M(t)`$ approach zero as the system evolves in time, as expected from winding mode annihilation (see Fig. III). As a result loop production continues until all the winding modes are annihilated at which time $`g(t)`$ approaches a constant (see Fig. III). Note from Figure III that before converging to zero, the two winding numbers converge to eachother. This is a necessary condition for isotropization. This result seems reasonable, since the rate of annihilation depends on the number of winding states present. Thus, as the number of winding modes decreases in one direction, the other direction should have an larger rate of annihilation. This suggests a sort of equilibration process due to the presence of winding modes.

However, the above alone is not enough evidence to prove that isotropization of the dimensions occurs. To prove that isotropization occurs we need to study the geometry of the background. We can see from Figure III that the expansion rates, $`l(t)`$ and $`q(t)`$, converge at late times. Note, however, that at the same time the ratio of scale factors continues to increase. A more quantitative definition of isotropization can be obtained by defining the average Hubble parameter $`\overline{\lambda }`$ and the anisotropy parameter $`𝒜`$ ,

$`\overline{\lambda }{\displaystyle \frac{1}{D}}{\displaystyle \underset{i=1}{\overset{D}{}}}\lambda _i`$ (38)
$`𝒜{\displaystyle \frac{1}{D}}{\displaystyle \underset{i=1}{\overset{D}{}}}{\displaystyle \frac{(\lambda _i\overline{\lambda })^2}{\overline{\lambda }^2}}.`$ (39)
In our case we have $`D=3`$ and the anisotropy parameter becomes,

$$𝒜=2\frac{(l(t)q(t))^2}{(2l(t)+q(t))^2}.$$
(40)
We find that for any amount of initial anisotropy, $`𝒜`$ approaches a maximum value and then goes to zero, as can be seen from Figure III. Furthermore, for the case of isotropic initial expansion but inequivalent winding mode numbers we find that $`𝒜`$ reaches a larger maximum but the conclusion of isotropization at later times is unchanged.

## IV Conclusion

We have generalized the equations of BGC to the anisotropic case including the effects of winding state annihilation and loop production. We address the issue of isotropization of the three large dimensions quantitatively by introducing the anisotropy parameter $`𝒜`$. Our analysis indicates that for an arbitrary amount of initial anisotropy, the anisotropy will reach a maximum early in the evolution and then approach zero at later times. Thus, we have explained how isotropy arises as a natural consequence of BGC.

###### Acknowledgements.

RB was supported in part by the U.S. Department of Energy under Contract DE-FG02-91ER40688, TASK A. SW was supported in part by the NASA Graduate Student Research Program. RB also wishes to thank the CERN Theory Division and the Institut d’Astrophysique de Paris for hospitality and financial support during the time the work on this project was carried out.