text:
# The study of spin-spin correlations in quasi-one-dimensional Heisenberg antiferromagnetic clusters

## Abstract

For a class of quasi-one-dimensional clusters, by using exact diagonalization, we study the effect of side spins on the spin-spin correlations on chain. Our calculations show that the side spins added in the same sublattice can effectively strengthen the spin-spin correlations in large distance region and make the change tend to flat. It is exactly proved that periodically adding side spins can set up magnetic long-range orders in the ground state. Also we investigate the effect of the density of side spins on correlation strength. The case that two sublattices have different localized spins is discussed.

Spin-spin correlation and magnetic long-range order (LRO) are of fundamental importance for quantum spin systems. They have been studied by many analytical and numerical works. For the bipartite Heisenberg antiferromagnetic (AF) systems, the spin-spin correlations of any two sites in the ground state (GS) always are AF, i.e. the correlation functions are positive when two sites belong to the same sublattice whereas negative the different sublattices. However, it is not sufficient to set up the AF LRO. The dimensionality of system is one of the key factors. It has been rigorously proved that there exists Néel order in the GS  of the three-dimensional system, whereas no Néel order for the one-dimensional one. Although there is no magnetic LRO for one- and two-dimensional (1D and 2D) Heisenberg antiferromagnets at temperature $`T>0`$ by Mermin-Wager theorem , $`T=0`$ may be the critical point. When the localized spin $`S1`$, it has been proved that the GS of the 2D Heisenberg antiferromagnet has Néel order . Notwithstanding the rigorous proof of Néel order when $`S=\frac{1}{2}`$ has not been established until now, many numerical and analytical works support its existence at $`T=0`$ .

The situation of quasi-one-dimensional (QOD) systems may be diverse due to their various geometric structures. The discovery of a heavy-fermion phenomenon in the $`Ce`$-doped Neodymium cuperate  has led to an increasing interest in the study of strongly correlated electrons coupled antiferromagnetically to magnetic moments . W. Zhang et al. have investigated the case of a single magnetic impurity and found that the spin-spin correlation function between the impurity spin and spins in the chain extends over long range . Obviously, the magnetic LRO can not be established through adding a single impurity spin to an 1D chain. The situation may change when impurity spins are periodically added to the chain in a special way. In this paper, by using exact diagonalization, we investigate a class of QOD Heisenberg AF clusters, in which there are some impurity spins sit beside an 1D finite chain (we call them as side spins), and explore the effects of side spins on spin-spin correlations. Our numerical results indicate that the spin-spin correlations in the region of large distances can be enhanced by adding side spins in the same sublattice. When periodically adding side spins, the decay of spin-spin correlations with distances slows down and becomes obviously flat in the range of large distances. For infinite 1D chains, it is analytically proved that adding side spins can set up magnetic LROs. Also, we investigate the variation of magnetic LROs with the density of side spins, and find that the decay of ferromagnetic (F) LRO is faster than that of AF LRO.

We investigate the spin-$`\frac{1}{2}`$ Heisenberg AF system with interactions of nearest neighbors, whose Hamiltonian reads

$$H=J\underset{(𝐢,𝐣)}{}\stackrel{}{S}(𝐢)\stackrel{}{S}(𝐣)\text{,}$$
(1)
where $`J>0`$. $`(𝐢,𝐣)`$ denotes the sum over pairs of nearest neighbors. The spin-spin correlation can be written as $`\mathrm{\Delta }(R_iR_j)G\left|\stackrel{}{S}(𝐢)\stackrel{}{S}(𝐣)\right|G`$. Here, $`R_i`$ and $`R_j`$ are the coordinates of sites $`i`$ and $`j`$ respectively, and $`|G`$ represents the GS. Let $`\overline{i}`$ denote the site by the $`i`$th site of the chain. At first, we consider a finite chain of $`23`$ sites under free boundary conditions and study the effect of a single side spin on the spin-spin correlations. By exact diagonalization, we calculate $`\mathrm{\Delta }(R_iR_j)`$ for site 2 (i.e. $`j=2`$ and $`i=3,4,5,\mathrm{}`$), and the numerical results are shown in Fig. 2, here the side spin is placed at site $`\overline{2}`$, $`\overline{3}`$ and $`\overline{4}`$ respectively (see Fig 1.(a)). When the side spin is added on site $`\overline{3}`$, $`\mathrm{\Delta }(R_iR_j)`$ becomes weaker than that of the pure 1D chain. Namely, the side spin weakens the spin-spin correlations. But, when the side spin is at site $`\overline{2}`$ (or $`\overline{4}`$), $`\mathrm{\Delta }(R_iR_j)`$ becomes stronger in the region of large distances. In this case, the side spin strengthens the spin-spin correlations for most of $`(R_iR_j)`$. When a single side spin is placed at sites $`\overline{5}`$, $`\overline{7}`$, $`\overline{9}`$ and $`\overline{11}`$ respectively, our calculations show that the effect of the side spin is similar to that of the side spin at site $`\overline{3}`$. And, when the side spin is placed at site $`\overline{6}`$, $`\overline{8}`$, $`\overline{10}`$ and $`\overline{12}`$ respectively, it is similar to that of the side spin at site $`\overline{2}`$.

Under the periodic boundary conditions, the case of a single side spin has been studied by Monte Carlo simulations and the spin correlation function is found to extend over long range . For such a system, the GS always takes the lowest possible total spin (LTS) no matter where the side spin is attached. But the situation for a finite chain under free boundary conditions becomes little complicated. The finite chain can be divided into two ‘sublattices’ as doing for infinite systems. Sites $`1`$, $`3`$, $`5`$, …, $`21`$ and $`23`$ belong to sublattice $`A`$, and sites $`2`$, $`4`$, $`6`$, …, $`20`$ and $`22`$ sublattice $`B`$ (Fig. 1(a)). Adding the single side spin at different sites can leads to the GS taking different values. When the side spin sits by the site of sublattices $`B`$, it belongs to sublattice $`A`$. By Lieb-Mattis theorem , the GS has the global spin 1 and is 3-fold degenerate. In other words, the GS takes a higher total spin than its lowest possible value. But, when the side spin sits by the site on sublattices $`A`$, it belongs to sublattice $`B`$. The GS takes the LTS zero. Our calculations show that the effect of a single side spin depends on where the side spin is. In other words, the total spin of the GS is related to the behaviors of spin-spin correlations. One important question concerned is how the spin correlations change when more side spins are added. When the side spins added are in the same sublattice, the total spin of the GS becomes higher, and the difference between it and the lowest possible total spin larger. From the above calculations, we speculate that $`\mathrm{\Delta }(R_iR_j)`$ will become stronger in the region of large distances.

For further details, we consider the cases of two side spins. If adding two side spins at site $`\overline{2}`$ and $`\overline{4}`$, the GS has the global spin $`\frac{3}{2}`$ while the lowest possible total spin is $`\frac{1}{2}`$. But, if adding two side spins by site $`2`$ and $`3`$, the GS takes the lowest possible total spin $`\frac{1}{2}`$. From the above notion, one can speculate that spin-spin correlations will become stronger if adding two side spin at sites $`\overline{2}`$ and $`\overline{4}`$, whereas those weaker at sites $`\overline{2}`$ and $`\overline{3}`$. We calculate the spin-spin correlations between site $`2`$ and others. The numerical results are shown in Fig. 2. They agree with the above speculation.

Also, we calculate spin-spin correlations for the following three cases (see Fig. 1(a)): 1). 4 side spins sitting at sites $`\overline{2}`$, $`\overline{8}`$, $`\overline{14}`$ and $`\overline{20}`$ on a finite chain of 21 sites; 2). 5 side spins sitting at sites $`\overline{2}`$, $`\overline{6}`$, $`\overline{10}`$, $`\overline{14}`$ and $`\overline{18}`$ on a finite chain of 19 sites; 3). 8 side spins sitting at sites $`\overline{2}`$, $`\overline{4}`$, $`\overline{6}`$, $`\overline{8}`$, $`\overline{10}`$, $`\overline{12}`$, $`\overline{14}`$ and $`\overline{16}`$ on a finite chain of 17 sites. Their ground states have the global spins $`S^T=\frac{5}{2}`$, $`3`$ and $`\frac{9}{2}`$ respectively. The numerical results are plotted in Fig. 3. Our calculations show that periodically adding side spins in the same sublattice can obviously enhance the spin-spin correlations in the large distance region. And, as the density of side spins increasing, spin-spin correlations become stronger and decay more slowly. Approximately, we can fit the spin-spin correlations in exponential way, i.e. $`\mathrm{\Delta }(R_iR_j)=C\mathrm{exp}(\left|R_iR_j\right|/\zeta )`$, here $`C`$ is coefficient and $`\zeta `$ the correlation length. We calculate $`\zeta `$ by $`\mathrm{\Delta }(R_9R_2)`$ and $`\mathrm{\Delta }(R_{15}R_2)`$. For the bare chain of 21 sites, $`\zeta =9.30`$. But, for the above three cases, $`\zeta =51.38`$, $`77.27`$ and $`7345.38`$. It shows that the large distance behavior of spin-spin correlations is enhanced by adding side spins. Especially, $`\zeta 10^3`$ for the third case. It is much larger than the correlation length of bare chain, even the size of system. Consequently, side spins can slow down effectively the decay and make the variation become flat in the region of large distances (Fig. 3). It seems to exist the AF LRO.

Now, we turn to consider an infinite QOD system, which is constructed by periodically adding side spins in the same sublattice (Fig. 1(b)). The total number of sites is $`N=K(l+2)`$, here $`K`$ denotes the number of cells and $`l`$ the number of sites on chain between every two side spins. $`l`$ must take odd, i.e. $`l=2k+1`$, here $`k=0,1,2\mathrm{}`$. Supposing the side spin is in sublattice $`A`$, the number of sites of sublattice $`A`$ is $`N_A=K(l+3)/2`$ and that of sublattice $`B`$ $`N_B=K(l+1)/2`$. By the Lieb-Mattis theorem , the global spin of the GS is

$$\mathrm{\Lambda }=K\left|(k+2)S_A(k+1)S_B\right|\text{,}$$
(2)
where $`S_A`$ and $`S_B`$ are the values of localized spins on sublattice $`A`$ and $`B`$ respectively. To investigate the existence of magnetic LRO, one needs to calculate the quantity

$$g(𝐪)G\left|\stackrel{}{S}(𝐪)\stackrel{}{S}(𝐪)\right|G\text{.}$$
(3)
$`\stackrel{}{S}(𝐪)=\frac{1}{\sqrt{N}}_𝐣\stackrel{}{S}(𝐣)\mathrm{exp}(i𝐪𝐣)`$, $`𝐪`$ is a reciprocal vector. The criterion of magnetic LROs is that $`g(𝐪)O(N)`$ at some $`𝐪`$. If $`g(𝐐)O(N)`$, here $`𝐐=(\pi ,\pi ,\mathrm{},\pi )`$, there exists AF LRO. If $`g(\mathrm{𝟎})O(N)`$, there exists F LRO. Following the approaches developed by G. S. Tian , one can obtain

$$g(𝐐)>g(\mathrm{𝟎})=\frac{\mathrm{\Lambda }^2+\mathrm{\Lambda }}{N}.$$
(4)
We discuss the case that the sites on two sublattices have the equal localized spin, i.e. $`S_A=S_A=S`$. From equation (2), we readily obtain $`\mathrm{\Lambda }=NS/(l+2)`$. As long as $`l`$ is finite, it is always true that

$$g(𝐐)>g(\mathrm{𝟎})>O(N).$$
(5)
From these inequalities, one can conclude that AF and F LROs coexist in the GS, and the former is predominant. In other words, side spins set up magnetic LROs.

Obviously, although there always exist magnetic LROs for finite $`l`$, both of F and AF correlation strengths will depend on the density of side spins, which is defined as $`\eta 1/(l+2)`$. We introduce the following two quantities to measure F and AF correlation strengths respectively,

$$\mathrm{\Gamma }_F\frac{1}{N^2}\underset{𝐢,𝐣}{}G\left|\stackrel{}{S}(𝐢)\stackrel{}{S}(𝐣)\right|G$$
(6)
and

$$\mathrm{\Gamma }_{AF}\frac{1}{N^2}\underset{𝐢,𝐣}{}\lambda _{\mathrm{𝐢𝐣}}G\left|\stackrel{}{S}(𝐢)\stackrel{}{S}(𝐣)\right|G,$$
(7)
where $`\lambda _{\mathrm{𝐢𝐣}}=1`$ when sites $`𝐢`$ and $`𝐣`$ belong to the same sublattice and $`\lambda _{\mathrm{𝐢𝐣}}=1`$ when sites $`𝐢`$ and $`𝐣`$ the different sublattices. $`\mathrm{\Gamma }_F`$ and $`\mathrm{\Gamma }_{AF}`$ will decrease as $`l`$ increasing. When $`l`$ reaches to the infinite, the system changes into an 1D chain and the magnetic LROs vanish (i.e. $`\mathrm{\Gamma }_F=\mathrm{\Gamma }_{AF}=0`$). Possibly, the variation speeds of $`\mathrm{\Gamma }_F`$ and $`\mathrm{\Gamma }_{AF}`$ are different. We calculate the ratio $`\rho \mathrm{\Gamma }_F/\mathrm{\Gamma }_{AF}`$ for $`l=1`$, $`3`$ and $`5`$, and plot the data in Fig. 4. $`\rho `$ decreases as $`\eta `$ decreasing approximately in the linear way. Then the decay speed of F correlation strength is faster than that of AF correlation strength. In other words, the ferrimagnetism becomes weaker and weaker as the density of side spins decreasing.

Ref.  has proved that the GS of 1D Heisenberg AF chain has magnetic LROs when its two sublattices have unequal localized spins. It is interesting to investigate the QOD system with unequal localized spins. From equation (2) and (4), we can conclude that if

$$\frac{S_A}{S_B}=\frac{k+1}{k+2}\text{,}$$
(8)
$`g(\mathrm{𝟎})=0`$ due to $`\mathrm{\Lambda }=0`$. It means that there is no F LRO. The simplest case is $`S_A=\frac{1}{2}`$ , $`S_B=1`$ and $`l=1`$. One can give a spin picture of the GS in valence-bond version. The spin on sublattice $`B`$ can be divided into two $`\frac{1}{2}`$ spins. One of them forms a singlet with the nearest side spin, and the other combines with its nearest neighbor on chain. We think this kind of configurations governs the physics of the GS. And it is responsible for the F LRO to vanish. Although we have not exactly proved that the AF LRO can not exist in the GS of such systems, we believe it is true. But, if $`S_A/S_B(k+1)/(k+2)`$, the AF and F LROs coexist in the GS.

We would like to thanks Prof. G. S. Tian for suggestions. The part of calculations in this work have been done on the SGI Origin 2000 in the Group of Computational Condensed Matter Physics, National Laboratory of Solid State Microstructures, Nanjing University. This work was partially supported by the Ministry of Science and Technology of China under Grant No. NKBRSF-G19990646.

Figure Captions:

Fig. 1: (a) the 1D Heisenberg AF chain. (b) the QOD Heisenberg AF chain. There are $`l`$ sites between site $`P_1`$ and $`P_2`$, and $`l=2k+1`$, here $`k=0,1,2\mathrm{}`$.

Fig. 2: Spin-spin correlations between site 2 and others (see Fig. 1(a)) for these cases: no side spin (square); single side spin at site $`\overline{2}`$ (cross $`+`$), $`\overline{3}`$ (circle) and $`\overline{4}`$ (cross $`\times `$) respectively; two side spins at sites $`\overline{2}`$ and $`\overline{4}`$ (triangle); two side spins at sites $`\overline{2}`$ and $`\overline{3}`$ (diamond). Here, $`R=R_iR_2`$. The length of chain is 23.

Fig. 3: Spin-spin correlations between site 2 and others (see Fig. 1(a)) for these cases: bare chain of 21 sites (square); 4 side spins at sites $`\overline{2}`$, $`\overline{8}`$, $`\overline{14}`$ and $`\overline{20}`$ on a chain of 21 sites (diamond); 5 side spins by site $`\overline{2}`$, $`\overline{6}`$, $`\overline{10}`$, $`\overline{14}`$ and $`\overline{18}`$ on a chain of 19 sites (triangle); 8 side spins by site $`\overline{2}`$, $`\overline{4}`$, $`\overline{6}`$, $`\overline{8}`$, $`\overline{10}`$, $`\overline{12}`$, $`\overline{14}`$ and $`\overline{16}`$ on a chain of 17 sites (solid circle). Here, $`R=R_iR_2`$.

Fig. 4: $`\rho `$ vs. $`\eta `$. The solid line is obtained by fitting $`\rho `$ in linear way.


format: text

text:
# Breakdown of Energy Equipartition in a 2D Binary Vibrated Granular Gas

(November 20, 2001)

## Abstract

We report experiments on the equipartition of kinetic energy between grains made of two different materials in a mixture of grains vibrated in 2 dimensions. In general, the two types of grains do not attain the same granular temperature, $`T_g=1/2m<𝐯^\mathrm{𝟐}>`$. However, the ratio of the two temperatures is constant in the bulk of the system and independent of the vibration velocity. The ratio depends strongly on the ratio of mass densities of the grains, but is not sensitive to the inelasticity of grains. Also, this ratio is insensitive to compositional variables of the mixture such as the number fraction of each component and the total number density. We conclude that a single granular temperature, as traditionally defined, does not characterize a multi-component mixture.

When two different gases are put in thermal contact with each other they eventually reach thermal equilibrium, with equipartition of kinetic energy between the gas molecules and a single temperature for the system as mandated by the zeroth law of thermodynamics. In this article we present an experiment designed to test the validity of energy equipartition in a mixture of two kinds of macroscopic grains driven into a non-equilibrium steady state by external vibration. Kinetic theories for granular systems have long used the idea of a granular temperature defined as the average kinetic energy of the particles, $`T_g=1/2m<𝐯^\mathrm{𝟐}>`$; the value of $`T_g`$ in a given system is determined by a balance between the power input by external driving sources and the energy lost in interparticle collisions. The results of such calculations of transport in single-component systems have been validated in both simulations HelalBiben97 ; BizonShattuck99-1  and experiments RouyerMenon01RerichaBizon01YangHuan01 . It is natural to ask whether the variable $`T_g`$ obeys a zeroth law, and in particular, whether a mixture of grains is characterized by a single, shared, granular temperature.

The earliest theoretical developments of kinetic theories for binary mixtures of spheres FarrellLun86 ; JenkinsMancini89  or of discs JenkinsMancini87 , start with a clearly identified assumption that a single temperature variable, $`T_g`$, characterizes the entire mixture. Most subsequent work preserves this assumption. An important exception is the work of Garzo and Dufty GarzoDufty99  who study theoretically the “cooling” of a binary granular mixture from an initial distribution of particle velocities. They conclude that the cooling rate of each species in the mixture is the same, although their temperatures are different throughout the cooling process. Experiments and simulations on binary mixtures of granular materials have primarily focused on segregation and mixing of different types of grains OttinoKhakhar00 , and for the most part, have not considered the issue of equipartition. However, an experiment on a vibrated monolayer of grains by Losert et al. LosertCooper99  found that grains of two different sizes and masses do not attain the same value of $`T_g`$. In this article, we show that such differences persist in the bulk of a 2D system and that in general, the two components of a mixture arrive at a ratio of granular temperatures, $`\gamma 1`$. We explore the dependence of the temperature ratio, $`\gamma `$, on particle properties as well as on compositional parameters of the mixture. Our data argue for a redefinition of the granular temperature in order to accommodate extension to a multi-component system.

We make binary mixtures of spherical glass balls with aluminum, steel or brass balls. The two relevant material properties being varied are mass and inelasticity (see Table 1), and the two compositional parameters being varied are the number fraction of each component, $`x`$, and the number density, specified in terms of the the average occupied area fraction, $`\rho _{avg}`$, of the particles in the cell. To determine the importance of particle material properties we hold average number density and number fraction fixed and compare mixtures of steel (heavy-elastic) and glass (light-elastic), brass (heavy-inelastic) and glass (light-elastic), and aluminum (light-inelastic) and glass (light-elastic) in terms of the $`T_g`$ of each component and the temperature ratio $`\gamma `$. Next, to determine whether failure of equipartition is a bulk effect, we hold the number fraction fixed for a steel-glass mixture and measure the temperature ratio as the number density of balls is varied. Finally, to test the dependence of the temperature ratio on the relative fraction concentration, we hold the number density fixed for a brass-glass pair, and measured the temperature ratio of the components as the number fraction of brass is varied.

The balls (diameter: $`d=1.600\pm 0.002mm`$) are confined to a vertical, rectangular, Delrin cage (32 $`d`$ high x 48 $`d`$ wide x 1.1 $`d`$ thick) sandwiched between two parallel plates of anti-static coated Plexiglas. An electromechanical shaker (LDS 500L) vibrates the cage vertically at a frequency of $`60Hz`$ and amplitudes up to $`2.4d`$, producing maximum accelerations, $`\mathrm{\Gamma }`$, and velocities, $`v_0`$, of $`56g`$ and $`1.45m/s`$ respectively. The motions of the balls are recorded with a monochrome high-speed camera (Kodak Motioncorder) at a rate of 2000 frames/s. Each ball is located with a precision of 0.03 $`d`$ ($`40\mu m`$). Glass balls are unambiguously distinguished from their metallic partners by using a combination of reflected and transmitted light, and by thresholding appropriately. Most of the results we discuss here are taken in a rectangular (10 $`d`$ x 20 $`d`$) window centered in the geometrical center of the cell.

The data presented in this article are all taken at relatively large values of excitation, corresponding to $`\mathrm{\Gamma }>32g`$ and $`v_0>0.86m/s`$. Beyond this scale of excitation, the vertical density profile does not vary much and tends to an asymptotic profile that is symmetric about the center of the cell, as shown in Fig. 1. The density profile is identical for two different species when their masses are matched, even if they have different inelasticity, as in the case of aluminum-glass mixture shown in Fig. 1(a). In contrast, the profiles of the other two mixtures (steel-glass and brass-glass) in Fig. 1(b), show a slightly higher concentration of heavy particles in the center of the cell and a more uniform distribution of light particles throughout the cell. However, in all three mixtures, while there are differences in the overall density distribution, the particles in the mixture are locally well-mixed, with no obvious tendency for balls of a particular type to cluster 222By contrast, at low accelerations, there is a clear tendency toward complete segregation, as in the simulations of S. McNamara and S. Luding, IUTAM symposium on Segregation in Granular Flows, eds. A. D. Rosato and D. L. Blackmore, (Kluwer Academic Publishers, Dordrecht, 2000), p. 305; S. Luding et al., ibid. p. 297. when vigorously excited as in this case. Finally, we have not observed any gradients of velocity or density in the horizontal direction.

The principal observation of our experiment is that the two components of a binary mixture do not always equilibrate to the same granular temperature. Figure 2(a) shows $`T_g`$ of three pairs of materials as a function of the squared vibration velocity of the cell, $`v_0^2`$, with the number density being held fixed at $`\rho _{avg}=0.096`$ and the number fraction of glass balls fixed at $`x=1/2`$. The granular temperature $`T_g`$ is proportional to the squared vibration velocity $`v_0^2`$ with the proportionality constant being a function of the inelasticity of the balls and number density $`\rho _{avg}`$. As shown in Fig. 2(b), $`T_g`$ equilibrates for the alum-glass mixture, but does not equilibrate for the steel-glass and brass-glass mixtures. In terms of particle properties, equipartition is achieved when the mass ratio is close to one ($`m_{glass}/m_{alum}=0.92`$), and fails when the masses are mismatched ($`m_{glass}/m_{steel}=0.29`$, $`m_{glass}/m_{brass}=0.33`$). Even when the temperatures of the two components do not equilibrate, the ratio of the temperatures do not vary as a function of $`v_0^2`$ within the error bars, as shown in Fig. 2(b). Furthermore, the ratio of temperatures are very close to one another ($`\gamma _{gs}=T_{glass}/T_{steel}=0.66\pm 0.06`$ and $`\gamma _{gb}=T_{glass}/T_{brass}=0.69\pm 0.09`$), even though the inelasticity of the metal partner in each case is quite different (Table 1). We conclude that the temperature ratio differs from one when the mass ratio departs significantly from one, but does not depend sensitively on the relative value of the inelasticity of the two components.

The failure of equipartition is not merely a boundary effect, nor is the value of $`\gamma `$ that we report in Fig. 2(b), an inhomogeneous average over different regions of the cell. This is shown in Fig. 3 where the temperature profiles of steel and glass are plotted as a function of the height of the cell, z\[$`d`$\], for $`v_0=1.45m/s`$. In a collision with a moving wall which is effectively infinitely massive, both species pick up the same velocity, which implies that the heavier species (steel) picks up greater kinetic energy. The temperature ratio $`\gamma `$ is thus different from one at the walls; this is analogous to the previous observations of Losert et al LosertCooper99 . Remarkably, however, the temperature ratio quickly reaches a constant value ($`\gamma =0.69`$) in the interior of the cell, even as the $`T_g`$ of each component continues to decrease toward the middle of the cell as a result of losses from inelastic collisions. Therefore the effect of inelastic collisions is to drive the value of $`\gamma `$ in the interior of the cell toward a value different from one, which as we mentioned above, depends on the mass ratio of components, but is insensitive to their inelasticity, and to the value of $`v_0`$.

We pursue further the possibility that the value of $`\gamma `$ obtained is due to boundary effects by increasing the number density, $`\rho _{avg}`$, of a steel-glass mixture holding the number fraction fixed at $`x=1/2`$. While this is not exactly equivalent to increasing the system size, it does tune the ratio of system size to mean free path, so that we vary the number of collisions experienced by balls in transferring energy from the boundary to the interior. The temperature ratio, $`\gamma =T_{glass}/T_{steel}`$ is plotted against of the square vibration velocity, $`v_0^2`$, in Fig. 4, for mixtures whose number density, $`\rho _{avg}`$, varies by a factor of 4. Once again, the ratio of the two temperatures are essentially independent of $`v_0^2`$. Importantly, the temperature ratio maintains an average value of $`\gamma =0.60\pm 0.06`$, with no obvious systematic trend as $`\rho _{avg}`$ changes. This result further confirms that the failure of equipartition is a bulk effect, and also suggests that the value of $`\gamma `$ is insensitive to the number density of the mixture.

Finally, we study the effect on equipartition of varying number fraction in a mixture at fixed density. Figure 5 shows the temperature ratio, $`\gamma =T_{glass}/T_{brass}`$, of a brass-glass mixture as a function of $`v_0^2`$, for three different number fractions of brass, $`x`$: 1/2, 1/4, and 1/8. The temperature ratios for the different relative concentrations are all the same within experimental errors. This demonstrates that the components do not reach thermal equilibrium, even though the relative number of collisions with glass balls increase for a given brass ball, as the number fraction of brass is decreased. It has been predicted SantosDufty01  that equipartition does not generally occur even in the tracer limit of a single foreign particle in a sea of other particles. This is consistent with our observation that the temperature ratio $`\gamma `$ does not depend on number fraction.

As mentioned earlier, Garzo and Dufty GarzoDufty99  study a freely-cooling binary system, and find a temperature ratio dependence on material parameters as well as on compositional parameters. Since there are no comparable predictions for a driven system, we have compared our data with these predictions and find quite good agreement with our results as a function of material parameters. However, they predict a strong dependence on number fraction of the components, where none is observed. The comparison between theory and experiment encourages further theoretical development in driven systems.

In conclusion, we have observed that energy equipartition does not generally hold for a binary vibrated granular system. However, the ratio between the steady state values of the granular temperatures of the two components does robustly satisfy some general trends. The temperature ratio is constant in the interior of a system, and is insensitive to vibration velocity and inelasticity, while depending strongly on mass density. Furthermore, the ratio is insensitive to compositional parameters such as number fraction of each component and average number density. This opens the possibility that just as in the work of Garzo and Dufty GarzoDufty99  for a cooling system, there is a well-defined and calculable ratio of temperatures achieved in a steady state system. Thus a modified definition of granular temperature with a prefactor depending chiefly on material properties, may adequately describe a multi-component system. Numerical simulations could play a useful role, in that material properties may be continuously tuned, as opposed to experiments in which only a limited pool of materials is readily available.

###### Acknowledgements.

We gratefully acknowledge support from NSF DMR CAREER 987433, and ACS PRF G33730. We also thank J. L. Machta for useful conversations and F. Rouyer for contributions to this work.


text:
# Strangeness Production as a Diagnostic Tool for Understanding Heavy Ion Reactions

## 1 Introduction

Hadronization of parton jets and multi-partonic systems is not directly calculable using Quantum Chromodynamics (QCD). This non-perturbative physics can only currently be described with phenomenological models that have been developed over the last 30 years. There are many such models including string fragmentation models, quark coalescence models, and statistical models. For the purposes of this presentation, we will use the statistical model as a baseline for comparing different colliding systems in terms of strangeness production. This is not meant to advocate one schematic description of hadronization over another.

Calculations assuming thermal and chemical equilibrium in hadron gas fireballs assume a universal hadronization mechanism and find a parton-hadron transition at a temperature of $`T170`$ MeV. The agreement between experimentally observed hadron yields in electron-positron jet fragmentation and in proton-(anti)proton reactions is reasonable. However, it is observed that particles carrying strange quarks are not in complete chemical equilibrium. There appears to be an additional suppression as expressed below:

$$\lambda _s=\frac{s\overline{s}}{\frac{1}{2}(u\overline{u}+d\overline{d})}$$
(1)
There appears to be a deficit of strange hadrons that is common in inclusive proton-proton, proton-antiproton, $`e^+e^{}`$ collisions as shown in Figure 1.

It was this observed suppression that led to the postulation that one could bring the strange hadrons into chemical equilibrium via pre-hadronization partonic interactions or post-hadronization re-scattering. However, the time scale for total strangeness ($`\mathrm{\Lambda }`$, $`K`$) equilibration via hadronic re-scattering is 10-100 fm/c, while for strange antibaryons ($`\overline{\mathrm{\Lambda }},\overline{\mathrm{\Xi }},\overline{\mathrm{\Omega }}`$) take more than 1000 fm/c. It is interesting to note that some recent calculations differ substantially from the above rate calculations.

The strangeness equilibration time in a deconfined system of a quark-gluon plasma is substantially shorter, of order 3-6 fm/c. This is due to contributions from gluon-gluon fusion to $`s\overline{s}`$ pairs. This time is well within the expected heavy ion collision lifetime of 10-15 fm/c before particles are free streaming. It is comparable to the time in heavy ion collisions that one might naively expect a deconfined region to exist before it cools to the point of hadronization.

## 2 Proton-Proton

Often the assumption is made that there is no quark-gluon plasma formed in proton induced reactions. This raises the question of the definition of the quark-gluon plasma. One answer is that it is an extended region of space with strong color fields characterized by a large density of deconfined partons and restoration of approximate chiral symmetry. Although the parton density achieved in elementary particle collisions is quite large, the volume is not, and one can argue whether it represents a well-defined thermodynamic state.

In fact, Bjorken speculated that in the “interiors of large fireballs produced in very high-energy $`p\overline{p}`$ collisions, vacuum states of the strong interaction are produced with anomalous chiral order parameters.” He applied a schematic picture referred to as the Baked Alaska in describing $`p\overline{p}`$ collisions, and the MINIMAX experiment at Fermilab searched for disoriented chiral condensates. Although MINIMAX yielded a null result, this picture is not ruled out.

Experiment E735 at the Fermilab Tevatron studied $`p\overline{p}`$ collisions to look for quark-gluon plasma signatures. Both E735 and UA1 and others observed a substantially larger source volume in high multiplicity $`pp`$ ($`p\overline{p}`$) events via two particle correlations measurements. They also observed transverse momentum spectra that have a mass ordering that some described via hydrodynamic flow. E735, UA1 and the AFS experiment at the ISR observed a substantial enhancement of strange particle production as a function of the total particle multiplicity. These observations led some to speculate that a quark-gluon plasma was being formed. However, these effects can also be explained consistently by a picture of event bias from the multiplicity selection on events with hard processes, in particular gluon jets. Thus, one needs to take care in comparing non-inclusive $`pp`$ data.

It is also interesting that statistical models often use the grand canonical ensemble. One can use the GCE even when the energy and other quantum numbers are conserved. The temperature and chemical potentials simply reflect the characteristics of the system. However one must be careful in the interpretation in particular of fluctuations quantities. If the volume of the system is large GCE is appropriate; however, for small volumes you must conserve quantum numbers (for example strangeness) in every event. Thus the canonical ensemble is relevant. In the CE, strangeness is suppressed in very small volume reactions and reaches the GCE limit only for large volume. It would be interesting to see if this can be tested in $`pp`$ events with different HBT radii. There may be technical difficulties in dealing with auto-correlations with HBT radii and other collision characteristics.

## 3 Proton-Nucleus

Inclusive $`pA`$ reactions have only a few binary $`NN`$ collisions, and are thus not so different than $`pp`$. However, measuring recoil nucleons in $`pA`$ reactions allows one to select events where the probability of many binary collisions is large. This type of analysis is now being done by experiment E910 at the BNL-AGS and NA49 at the CERN-SPS.

E910 observes that the $`K^+`$ and $`\mathrm{\Lambda }`$ production increases rapidly as a function of the number of collisions suffered by the proton and then appears to level off with no increase in production after $`4`$ collisions as shown in Figure 2.. The data do not appear to scale with the number of participating nucleons (as is often used as a baseline for comparison in AA collisions). E910 has pointed out that the difference between the $`\mathrm{\Lambda }`$ production data and participant scaling appears to increase linearly with the number of binary collisions up to three collisions. Three step valence quark stripping mechanisms and baryon junction models have been suggested to explain the data. The $`\mathrm{\Xi }`$ is also measured and shows an even stronger enhancement relative to participant scaling that the singly strange hadrons. Also, most of the additional production is at target rapidity. Thus, measurements only at mid-rapidity are missing most of the increased production. $`K^{}`$ production increases with $`N_{coll}`$, but after 4-5 collisions it sharply decreases. This decrease is potentially explained in that when the proton strikes the center of the nucleus ($`N_{coll}>5`$), there is a high probability for processes such as $`K^{}+N\pi +\mathrm{\Lambda }`$ to effectively absorb the $`K^{}`$.

NA49 also observes a similar enhancement for all strange particles, in particular the the doubly strange baryons. NA49 and WA97 observe that the $`\mathrm{\Xi }`$ yields increase faster than would be expected from simple participant scaling. Understanding the physical mechanism for the strangeness enhancement at both AGS and SPS energies is a theoretical challenge.

## 4 Nucleus-Nucleus

Heavy ion collisions at the AGS show an increase in the $`K^+/\pi ^+`$ ratio as a function of centrality or the number of participating nucleons, and thus total strangeness is enhanced. Most hadron yields can be described by a statistical equilibrium model within systematic errors, and one can see in Figure 1 the value for $`\lambda _s`$. At the lowest energies run at the AGS ($`p_{beam}=2GeV`$), strangeness production is near threshold and the yield drops drastically as expected.

Experiment E917 measures a ratio $`\frac{\overline{\mathrm{\Lambda }}}{\overline{p}}=3.6_{1.8}^{+4.7}`$, which confirms an earlier indirect measure by experiments E864 and E878. This large enhancement of strange antibaryons is above the statistical equilibrium level and may result from some balance of enhanced production of $`\overline{\mathrm{\Lambda }}`$ and annihilation losses for the $`\overline{p}`$. Unfortunately there is no measurement of $`\overline{\mathrm{\Xi }}`$. This is an excellent possible future measurement at the Japanese Hadron Facility (JHF) or an energy upgraded GSI.

Many experiments at the CERN-SPS also observe enhancement of strangeness in many channels. In particular, WA97 and NA49 observe a very large increase in the production of $`\mathrm{\Xi }`$, $`\mathrm{\Omega }`$, and $`\overline{\mathrm{\Xi }}`$, $`\overline{\mathrm{\Omega }}`$. The enhancement of participant scaling increases with increasing strangeness content $`E_\mathrm{\Lambda }<E_\mathrm{\Xi }<E_\mathrm{\Omega }`$. This is the same trend seen in $`pA`$. The statistical models give good agreement with a temperature of 170 MeV, in agreement with that obtained from high energy $`pp`$, $`p\overline{p}`$ and $`e^+e^{}`$.

The press release from CERN on the discovery of the quark-gluon plasma states that, “Since the hadron abundances appear to be frozen in at the point of hadron formation, this enhancement signals a new and faster strangeness-producing process before or during hadronization, involving intense re-scattering among quarks and gluons.” It seems that this conclusion is made without full consideration of data at lower energy and in $`pA`$ reactions.

At RHIC there is already preliminary data from all four experiments on strange and non-strange hadron yields. They also give a reasonable fit by a statistical model assuming equilibration of strangeness with T = 170 MeV. The enhancement of strangeness appears remarkably similar at the AGS, CERN-SPS and RHIC. I predict that LHC will see similar results.

## 5 Conclusions

Almost everywhere we look we observe enhancement of strange particle production. This effect is magnified as we go from proton-proton to proton-nucleus to nucleus-nucleus collisions. Many in the field assume this implies that strangeness production cannot be a signal of quark-gluon plasma formation. I believe this conclusion is premature. We must understand the global hadronization process from a region of vacuum with strong color fields in order to make progress. Future data at RHIC in $`pp`$, $`pA`$ and $`AA`$ are eagerly awaited.

## Acknowledgments

Useful discussions with Ken Barish, Sean Kelly, Dave Morrison and X.N. Wang are acknowledged.


text:
# A comment on the colour-colour diagrams of Low Mass X-ray Binaries

## 1 Introduction

Low mass X-ray binaries (LMXBs) hosting a neutron star can be observationally divided into two main categories, dubbed atolls and Z sources (Hasinger & van der Klis 1989). This classification is based on changes in both spectral and timing properties as the source varies. Z sources are named after a Z-shaped track they produce on an X-ray colour-colour diagram. Atolls can fall into one of the two spectral states, a hard, low luminosity ‘island’ or a soft, high luminosity ‘banana’. They trace a U-shaped or C-shaped track as the source spectrum evolves between the island and the banana (see e.g. fig. 1 in Méndez et al. 1999). These differences between the two LMXBs categories probably reflect differences in both mass accretion rate, $`\dot{M}`$, and magnetic field, $`B`$, with the Z sources having high luminosity (typically more than 50 per cent of the Eddington limit) and magnetic field ($`B10^9`$ G) while the atolls have lower luminosity (generally less than 10 per cent of Eddington) and low magnetic field ($`B10^9`$ G) (Hasinger & van der Klis 1989).

Both atolls and Z sources move along their tracks on the colour-colour diagram and do not jump between the track branches. Most of the X-ray spectral and timing properties, e.g. the kilohertz quasi-periodic oscillation (QPO) frequency (Méndez & van der Klis 1999), depend only on the position of a source in this diagram. This is usually parameterized by the curve length, $`S`$, along the track.

This strongly suggests that a single parameter determines the overall properties of the source, generally believed to be the accretion rate, which increases with $`S`$ from the horizontal (top) to the flaring (bottom) branch in the Z sources, and from the island to the banana in the atolls. The increasing characteristic frequencies in the power spectra along the track are then generally explained by the inner disc radius decreasing as a function of $`S`$ (e.g. the review by van der Klis 2000). However, the situation is more complex as $`S`$ is not simply related to the observed X-ray luminosity as would be expected if it is determined by the mass accretion rate (e.g. van der Klis 2000; 2001).

In this letter we present a compilation of RXTE data from three transient atoll sources. Their large amplitude of luminosity variation allows us to plot the full track in the colour-colour diagram, which appears to form a shape of ‘Z’, similar to Z sources.

## 2 The data

We have analysed the RXTE observations of the three atoll LMXBs: Aql X-1 (Cui et al. 1998; Reig et al. 2000), 4U 1608–52 and 4U 1705–44 (Hasinger & van der Klis 1989). 4U 1608–52 is a transient source and its most recent outburst in 1998 was observed by RXTE with a good coverage of both island and banana states. Aql X-1 is also a transient, showing outbursts in a time-scale of months to years. 4U 1705–44 is a strongly variable X-ray source, switching between the island and banana states on the time-scales of months.

We have used publicly available RXTE/PCA data of these three sources from PCA epochs 3 (between 1996 April 15 and 1999 March 22) and 4 (between 1999 March 22 and 2000 May 13). We selected data from detectors 0, 2 and 3, excluding all type I bursts and observations with very poor statistics. This gave 477 ks of data for Aql X-1, 219 ks for 4U 1608–52 and 172 ks for 4U 1705–44. The PCA light-curves for each energy channel were extracted in 128 s bins. These were used to build a colour-colour diagram, defining a soft colour as a ratio of 4–6.4 to 3–4 keV count rates, and a hard colour as a 9.7–16 over 6.4–9.7 keV ratio.

The response of the PCA detectors is slowly varying as a function of time. Additionally, the high-voltage settings of the instruments were altered between the PCA epochs 3 and 4. Therefore, the energy boundaries of each energy channel change in time, causing a shift in the colours. We have approximately taken these changes into account by reading these boundaries from the response matrices created for the beginning and end of each PCA epoch, and by linearly interpolating between them. When accumulating counts in each of the four energy bands (used for computing the colours), we have interpolated the number of counts for channels on the edges of these energy bands. We have checked the correctness of this procedure using Crab data from various observations in both PCA epochs. We computed colours and noticed that the position of Crab on the diagram still changes in time. This is due to the approximate character of the colours we had calculated. Therefore, to account for this variation, we have calculated multiplicative factors for the colours from the Crab data and applied them to our data (see e.g. Homan et al. 2001 for the similar method). The final result is presented in Fig. 1 (upper panels).

To make these plots clearer and to enhance the evolutionary track on the diagram we have rebinned the data using a nearest neighbour clustering technique. In each step of the iteration the two nearest bins (data points) on the colour-colour diagram were found. The total numbers of counts in both bins were added and a new count rate and colours were calculated. Thus, the two bins were replaced by one. The procedure was carried on until the assumed number of bins (between 8 and 10, as seen in the lower panels of Fig. 1) was reached. This method gathers together the data with similar spectral properties (colours) as opposed to increasing the length of the time bin, which averages data in time. The result is presented in Fig. 1 (lower panels). The size of the symbol is proportional to the logarithm of the count rate, so it gives an overview of the luminosity changes in the diagram.

## 3 Results

All three colour-colour diagrams in Fig. 1 show features characteristic of atoll sources that have been known for many years (e.g. van der Klis 1995; Méndez 1999). There is a banana in the lower part and several islands in the upper part of the diagram. The banana is significantly brighter than the islands, and the X-ray flux increases along the banana branch, from left to right.

However, this large compilation of data shows other features which have not previously been seen from sparser data sets. There are three distinct branches in the colour-colour diagram. The banana branch forms a lower horizontal track at hard colour of $``$ 0.4, while the island state begins along the diagonal track which connects to the left-hand end of the banana. But there is a part of the diagram not reported before: an upper horizontal branch, at hard colour of $``$ 0.7–0.8, connected to the upper right end of the diagonal branch. It is particularly pronounced in 4U 1705–44. This extends the previously known C-shaped pattern into a Z. There were hints of this extension in previous observations of 4U 1705–44 (Langmeier et al. 1989), but this is the first time that it is shown so noticeably.

To show this new branch more clearly we have created a combined colour-colour diagram of all three sources. Since these three atolls have similar spectral and timing properties, they probably have the same underlying accretion geometry and radiation mechanisms. Thus we might expect their spectral evolution, and hence the colour-colour diagrams to be the same, but despite similarities in shape, the colour-colour diagrams of Fig. 1 are shifted with respect to each other. This is mainly caused by differences in absorption, which strongly affects the soft colour, but there might also be subtle shifts in hard colour from differences in spin frequency (affecting disc-to-boundary layer luminosity ratio and hence the soft colour as well) and perhaps also the inclination angle. We have attempted to correct for these effects by simply shifting the diagrams in both colours. We have linearly transformed each 128 s data point of each source in a way that the left-hand edge of the lower branch have transformed colours of (1, 1) and the right-hand edge of the upper branch transforms to (2, 2). The resulting diagram of all three sources together is shown in Fig. 2. We can see that the three diagrams coincide very well, and the three branches forming a Z-shaped track are clearly visible.

We have analysed source movement along the Z-shaped track. For this, we have taken the original (i.e. not colour-binned) colour-colour diagram of 4U 1705–44 and traced how the position of the source (in 128 s data bins) moves with time. We confirm that motion in the diagram, including the upper branch, goes along the track and we do not notice jumps between the branches (though jumps cannot be completely excluded, since the data is rather sparse in time). 4U 1705–44 can cross the full width of the upper branch in about 10 days, while transition on the diagonal takes about five days (a similar transition in 4U 1608–52 took about three days). The movement in the lower branch is much faster, with time-scales of hours.

The lower panels in Fig. 1 show that the average X-ray flux at a given point on the track increases from left to right (i.e. with increasing soft colour) both in the upper and lower branches. To study the source flux along the track in details, for each grouped point in Fig. 1 we have extracted a typical PCA spectrum from a single RXTE pointing (one to three orbits, usually 3–10 ks of exposure) with the same colours and mean count rate. These spectra were fit with a model of a blackbody, Comptonized component and its reflection, with absorption set at 0.5, 1.5 and 1.2$`\times 10^{22}`$ cm-2, for Aql X-1, 4U 1608–52 and 4U 1705–44 respectively (Church & Bałucińska-Church 2001; Penninx et al. 1989; Predehl & Schmitt 1995). These show that the spectra all along the upper horizontal branch are rather hard, similar to those previously seen in the hardest island states. The bolometric flux was calculated by extrapolating the unabsorbed model spectrum, so it is somewhat model dependent but gives at least a zeroth order correction. In Fig. 3 we plot this bolometric flux in arbitrary units as a function of distance $`S`$ along the Z track. We have defined the right-hand edge of the upper branch as $`S=1`$, and the left-hand edge of the lower branch as $`S=2`$ and linearly interpolated to get the $`S`$ value of all the other data points (see Fig. 2). The bolometric fluxes are scaled so that they are roughly equal at the $`S=1`$ point. It is clear that the shapes of these curves are very similar for all 3 sources, showing a steady rise in the upper branch ($`S1`$), and then a drop on the diagonal branch, and then rising again along the lower branch.

This average bolometric flux most likely corresponds to the average mass accretion rate. Then, extending a LMXB paradigm we suggest that the accretion rate increases along the Z-shaped track in the colour-colour diagram. There is however a decrease of the flux on the diagonal branch. We speculate that this might associated with jet formation. It is known that jets are associated with state transitions in both neutron star and black hole transients (e.g. Fender & Kuulkers 2001). Alternatively, the mass accretion rate onto the central source (and hence the hard X-ray flux) could be reduced if much of the inflowing material was used to extend the disc inwards.

To convert X-ray fluxes into true luminosity requires a distance. For these three sources the best distance estimates are: 3.6 kpc for 4U 1608–52 from radius expansion of X-ray bursts (Nakamura et al. 1989), 4–6.5 kpc for Aql X-1 from optical spectroscopy of the companion star (Rutledge et al. 2001) and 6.3–8.2 kpc for 4U 1705–52 from modelling of X-ray bursts (Haberl & Titarchuk 1994). This gives the $`S=1`$ point at 5, 4–10 and 10–17 per cent of Eddington luminosity ($`L_{\mathrm{Edd}}=1.76\times 10^{38}`$ ergs s-1), for 4U 1608–52, Aql X-1 and 4U 1705-44, respectively. Since luminosity estimates from X-ray bursts are approximate and might contain unknown systematic errors (see e.g. in ’t Zand et al. 2001), it is possible that the position in the colour-colour diagram of all these three atolls depend on the X-ray luminosity in the same way, and that the $`S=1`$ point corresponds to $``$ 5–10 per cent of $`L_{\mathrm{Edd}}`$.

## 4 Discussion

The three atolls studied here are transient systems, where the luminosity changes by a factor of $`100`$, unlike the majority of atolls which are not very variable. The transients go down to very low luminosities where they show a new track on the colour-colour diagram which extends the previously known atoll (or C) shaped path into a Z. This new upper branch is clearly distinguished from a simple extension of the previously known C-shaped track. Below $`S=1`$, the track turns so that the average source luminosity decreases from right to left on the colour-colour diagram. We propose that all atolls would show such a track if their mass accretion rate could change by a large enough factor, but that their normal, rather small range in variability limits their observed colour-colour diagram to only a small section of the track.

A good example of this is 4U 1728–34, regarded as an archetypal atoll tracing a C-shaped track in the diagram (Méndez & van der Klis 1999; Di Salvo et al. 2001). However, the ratio of the highest to lowest count rate of this source is only about 2.5 (Di Salvo et al. 2001). Therefore, we suggest that the 4U 1728–34 data collected so far shows a colour-colour diagram which is limited to the diagonal and lower branches of a Z track. The source never goes to low enough luminosities to sample much of the upper branch.

Thus, we show that atolls share the same colour-colour topology with the Z sources. With increasing accretion rate they both trace out a similar Z-shaped pattern. Despite this similarity, these two LMXB categories are of course not the same. They differ in the luminosity range required to cover the whole Z track (Z sources are much less variable then the three transient atolls presented here) and in the time taken to trace out the upper and diagonal branches (Z sources move much faster). Another important difference is the luminosity at which the Z-shaped pattern arises. For the Z sources the transition from the upper (horizontal) branch to the diagonal (normal) branch occurs at around 0.5–1 $`L_{\mathrm{Edd}}`$, while for the atolls this transition occurs at luminosities about ten times less. There is also a significant difference in the spectral shape in the upper branch: spectra of the atolls are much harder, similar to hard state of black hole candidates.

The differences in luminosity and spectral shape can be reconciled in a model in which the fundamental difference between atolls and Z sources is magnetic field. Evolution along the Z track is caused by the increasing mass accretion rate, $`\dot{M}`$, decreasing the inner radius, $`R_{\mathrm{in}}`$, of a truncated disc. For the atolls the disc is probably truncated by evaporation (e.g. Różańska & Czerny 2000), leading to an inner, optically thin, hot flow which gives the hard X-ray spectrum. However, the evaporation efficiency decreases as a function of increasing mass accretion rate, so this cannot truncate the disc in the Z sources. Instead the truncation is likely to be caused by stronger magnetic field, but here the increased mass accretion rate means that the inner flow is much more optically thick, and so cooler.

We note that after this paper was submitted to MNRAS, another group independently presented very similar results (Muno, Remillard & Chakrabarty, astro-ph/0111370).

## Acknowledgements

We thank Didier Barret for attracting our attention to the 4U 1705–44 data. We also thank the referee, Michiel van der Klis, and Mariano Méndez for helpful discussions. This research has been supported in part by the Polish KBN grant 2P03D00514.


text:
# Constraining the Matter Power Spectrum Normalization using the SDSS/RASS and REFLEX Cluster surveys

## 1 Introduction

The present-day number density of galaxy clusters remains one of the most powerful constraints on the amplitude of matter perturbations in the Universe. This is usually defined in terms of $`\sigma _8`$, which is the dispersion of the mass field smoothed on a scale of $`8h^1\mathrm{Mpc}^1`$, where $`h`$ is the present value of the Hubble parameter in units of $`100\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$. Most often the local X-ray cluster temperature function is used for this purpose, given that the X-ray temperature has been the best observable from which to estimate cluster mass. While observation can readily give the X-ray temperature of clusters, theory can only easily predict the cluster mass function. To bridge the gap, theoretical modelling of clusters is used to provide a relation between mass and temperature, which in the most sophisticated treatments is taken to depend on both redshift and the underlying cosmology.

However, a drawback in the use of the local X-ray cluster temperature function is that only a few tens of clusters have had their X-ray temperature estimated. Consequently, authors using samples that only partially overlap have obtained significantly different cluster temperature functions, and thus estimates for $`\sigma _8`$ (e.g. compare Eke, Cole & Frenk 1996; Henry 1997, 2000; Markevitch 1998; Blanchard et al. 2000; Pierpaoli et al. 2001). Statistically it would be more robust if the number density of local clusters could be estimated from much larger samples.

We describe here a new approach which avoids working with the X-ray temperature function. We will instead use the luminosity function from the ROSAT-ESO Flux Limited X-Ray (REFLEX) galaxy cluster survey (Böhringer et al. 2001a,b), which contains 452 clusters, to estimate the local cluster number density. We relate X-ray luminosity to cluster mass by taking advantage of a weak shear lensing analysis (Sheldon et al. 2001) of a sample of 42 galaxy clusters identified in data from the commissioning phase of the Sloan Digital Sky Survey (SDSS), and cross-checked via correlation with the ROSAT All-Sky Survey (RASS) (Nichol et al. 2001). This relation is then used to compare the REFLEX luminosity function with the cluster mass function, leading to the estimation of $`\sigma _8`$, without any need to model the cluster mass–temperature relation (as in e.g. Viana & Liddle 1996, 1999).

## 2 Methodology

We consider spatially-flat cosmological models, where part of the energy density may be due to a cosmological constant, containing a spectrum of primordial adiabatic density perturbations. This family includes the current standard cosmological model, with a present-day matter density $`\mathrm{\Omega }_00.3`$, which provides the best fit to the full compilation of structure formation data (e.g. Durrer & Novosyadlyj 2001; Wang, Tegmark & Zaldarriaga 2001). We take the present-day shape of the matter power spectrum to be well approximated by that of a cold dark matter model with scale-invariant primordial density perturbations and shape parameter $`\mathrm{\Gamma }`$ in the range \[0.08,0.28\]; this interval is the average of the best-fit values coming from the preliminary analysis of the 2dF (Percival et al. 2001) and SDSS (Dodelson et al. 2001; Szalay et al. 2001) data, taking into account both statistical and systematic uncertainties. We will assume that $`\mathrm{\Gamma }`$ has an equal probability of taking any value within the interval given. Note, however, that ultimately the local cluster number density depends significantly only on $`\sigma _8`$ and $`\mathrm{\Omega }_0`$.

The expected halo mass function for each set of cosmological parameters is estimated via the fitting function of Jenkins et al. (2001), obtained by analysing data assembled from various large $`N`$-body simulations. Following White (2001) we consider the halo mass to be that given by the virial relation. The high-mass end of the halo mass function is better estimated using the Jenkins et al. (2001) result rather than the usual Press–Schechter ansatz (Press & Schechter 1974). The latter leads to a smaller number of high-mass halos at fixed $`\sigma _8`$, thus requiring systematically higher values of $`\sigma _8`$, by about eight percent, in order for the local cluster abundance to be reproduced. We cautiously model the uncertainty in the normalization of the mass function by means of a Gaussian distribution with a 10 percent dispersion around the mean (see Fig. 8 of Jenkins et al. 2001).

The local cluster number density was obtained by integrating the REFLEX X-ray luminosity function in the \[0.1,2.4\] keV band upwards from the mean luminosity of the 42 clusters in the SDSS/RASS sample, which is $`(0.17\pm 0.03)\times 10^{44}h^2\mathrm{erg}/\mathrm{s}`$. Their average redshift is 0.1, which roughly coincides with the mean redshift of the REFLEX clusters with similar or higher luminosity. For the 42 SDSS/RASS clusters, Sheldon et al. (2001) obtained, via statistical weak lensing analysis, a mean projected mass within $`r_{500}`$ of $`(0.9\pm 0.2)\times 10^{14}h^1M_{}`$, where $`r_{500}`$ is the radius within which the cluster mean density falls to 500 times the critical density at the redshift of observation. Following Sheldon et al. (2001) we assume the cluster density profile to behave like that of a singular isothermal sphere, i.e. $`\rho (r)r^2`$. We can now calculate the mean cluster mass within the three-dimensional radius $`r_{500}`$, and then convert this to a mean virial mass. Taking this mass as the lower limit in the integration of the Jenkins et al. (2001) mass function will then yield the expected local cluster number density for clusters more luminous than the mean of the SDSS/RASS sample, as a function of $`\sigma _8`$ and $`\mathrm{\Omega }_0`$. By comparison with the REFLEX estimate, best-fit values for $`\sigma _8`$ as a function of $`\mathrm{\Omega }_0`$ can then be obtained.

## 3 Results

Table 1 shows the 95 per cent confidence interval on $`\sigma _8`$ obtained using the REFLEX X-ray luminosity function and the relation between cluster mass and X-ray luminosity for the full SDSS/RASS sample. This interval was determined via Monte Carlo simulations, which incorporated all the uncertainties previously mentioned that affect the present estimation of $`\sigma _8`$. We find that the most probable value for $`\sigma _8`$ can be accurately represented by the fitting function

$$\sigma _8=0.38\mathrm{\Omega }_0^{0.48+0.27\mathrm{\Omega }_0},$$
(1)
with a 95 per cent uncertainty around 15%.

This result is significantly lower than, and barely compatible with, most determinations of $`\sigma _8`$ based on the local cluster X-ray temperature function (e.g. compare with Eke, Cole & Frenk 1996; Henry 1997, 2000; Viana & Liddle 1999; Blanchard et al. 2000; Pierpaoli et al. 2001), with the exception of Seljak (2001). His analysis differs from the others in that he used the relation between cluster temperature and mass derived by Finoguenov, Reiprich & Böhringer (2001) from X-ray data, rather than one obtained from hydrodynamical $`N`$-body simulations. The earlier work of Markevitch (1998) had already hinted at lower values for $`\sigma _8`$, at least in the case of $`\mathrm{\Omega }_00.3`$, if actual X-ray data was used to relate cluster temperature to mass. Similar results have been reached by Borgani et al. (2001), based on ROSAT Deep Cluster Survey (RDCS) data and the observed cluster X-ray temperature to luminosity relation, and Reiprich & Böhringer (2001), by means of an empirical cluster mass function derived using X-ray data from a large cluster sample.

Other methods of measuring $`\sigma _8`$ lead to conflicting results. While high-redshift Lyman-$`\alpha `$ forest analyses (Croft et al. 1999, 2000; McDonald et al. 2000), support our findings, estimates based on cosmic shear data (Hökstra, Yee & Gladders 2001; Maoli et al. 2001; Van Waerbeke et al. 2001) tend to favour higher values for $`\sigma _8`$ than those obtained here. Also, the 2dF galaxy survey, when combined with measurements of the amplitude of temperature anisotropies in the cosmic microwave background radiation (Efstathiou et al. 2001; Lahav et al. 2001), seems to require a value for $`\sigma _8`$ close to what we have found.

In order to test whether the low values obtained for $`\sigma _8`$ could be due to hidden systematic errors in the weak lensing method used for cluster mass estimation, we repeated the calculation of $`\sigma _8`$ for the two SDSS/RASS sub-samples discussed within Sheldon et al. (2001), namely

1. The 27 clusters with the lowest X-ray luminosities, on average $`(0.09\pm 0.02)\times 10^{44}h^2\mathrm{erg}/\mathrm{s}`$. They have a mean redshift of 0.09 and a mean projected mass within $`r_{500}`$ of $`M_{500}(1)=(0.7\pm 0.2)\times 10^{14}h^1M_{}`$.
2. The 15 clusters with the highest X-ray luminosities, on average $`(0.51\pm 0.04)\times 10^{44}h^2\mathrm{erg}/\mathrm{s}`$. They have a mean redshift of 0.17 and a mean projected mass within $`r_{500}`$ of $`M_{500}(2)=2.7_{1.1}^{+0.9}\times 10^{14}h^1M_{}`$.

Given that most of the REFLEX clusters with similar luminosities to those in the SDSS/RASS sample have a redshift between 0.05 and 0.2 (Böhringer et al. 2001a), with the higher luminosity clusters typically being at higher redshifts, we will assume that the REFLEX luminosity function provides a good representation of the underlying cluster luminosity function over this redshift interval. This is supported by an analysis of the Brightest Cluster Survey (BCS), which showed that there is no strong evidence for evolution in the cluster luminosity function out to at least $`z=0.2`$ (Ebeling et al. 1997).

Surprisingly, the most probable values for $`\sigma _8`$ according to each of these two sub-samples are rather different. Sub-sample 1 yields $`\sigma _8=0.37\mathrm{\Omega }_0^{0.52+0.32\mathrm{\Omega }_0}`$, almost indistinguishable from the result obtained from the full SDSS/RASS sample, while sub-sample 2 gives $`\sigma _8=0.50\mathrm{\Omega }_0^{0.47+0.24\mathrm{\Omega }_0}`$ which is substantially higher. This latter result is much more in line with standard $`\sigma _8`$ estimates from the local cluster X-ray temperature function.

## 4 Discussion

In this section, we discuss our results and present several tests to assess their robustness. First, to determine whether our results indicate an internal inconsistency in the weak lensing analysis of Sheldon et al. (2001), we performed Monte Carlo simulations where we used the data from sub-sample 1 to compute the inferred mean projected masses within $`r_{500}`$ for the clusters in sub-sample 2 and vice-versa, for various values of $`\mathrm{\Omega }_0`$. We included all the uncertainties previously mentioned in the simulations. We found that, provided $`\mathrm{\Omega }_0>0.3`$, sub-sample 1 implies a value for $`M_{500}`$(2) at least as large as that deduced by Sheldon et al. (2001) less than 5 per cent of the time, while sub-sample 2 implies a value for $`M_{500}`$(1) as small as that deduced by Sheldon et al. (2001) less than 10 per cent of the time. For values of $`\mathrm{\Omega }_0`$ between 0.1 and 0.3, the discrepancy decreases with $`\mathrm{\Omega }_0`$, but not significantly. We therefore conclude that, within the context of the cosmological models we discuss, the mean projected masses for the two SDSS/RASS cluster sub-samples presented in Sheldon et al. (2001) are only barely compatible within the uncertainties associated with the estimation of $`\sigma _8`$.

This inconsistency could result from the assumption that the cluster mass profile is that of a singular isothermal sphere (SIS), which was motivated by the fact that within the radius for which there was weak lensing data, the shear profile was close to that of a projected SIS. Although observationally the issue is still unresolved (e.g. Irwin, Bregman & Evrard 1999; Irwin & Bregman 2000; White 2000; De Grandi & Molendi 2001; Komatsu & Seljak 2001), it has become clear that $`N`$-body simulations produce clusters which on average have outer mass profiles, up to the cluster virial radius, that behave more like $`\rho (r)r^{2.5}`$ (Navarro, Frenk & White 1995, 1996, 1997; Thomas et al. 1998, 2001; Tittley & Couchman 2001). Consequently, we investigated the effect on the virial mass to X-ray luminosity relation of assuming an outer cluster mass profile parametrized by $`\rho (r)r^\beta `$. We varied $`\beta `$ between 2 (SIS case) and 3, and, as before, we took the mean projected masses within $`r_{500}`$ to be those given by Sheldon et al. (2001). This translates to fixing the inner projected cluster masses to be $`M_{500}`$(1) within $`0.39h^1\mathrm{Mpc}`$ for sub-sample 1, $`M_{500}`$(2) within $`0.57h^1\mathrm{Mpc}`$ for sub-sample 2, and $`0.9\pm 0.2\times 10^{14}M_{}`$ within $`0.41h^1\mathrm{Mpc}`$ for the full SDSS/RASS sample (Sheldon et al. 2001). As expected, we find that as $`\beta `$ increases, the cluster virial mass decreases, hence lowering $`\sigma _8`$. This change is greater for smaller $`\mathrm{\Omega }_0`$. Consequently, we find that the value for $`\sigma _8`$ derived using the Sheldon et al. (2001) and REFLEX data is barely affected by assuming different outer mass profiles if $`\mathrm{\Omega }_0`$ is close to one, whereas for low values of $`\mathrm{\Omega }_0`$ assuming a outer mass profile different from the true one may introduce a significant systematic error in the calculation of $`\sigma _8`$. For example, for $`\beta =3`$ we obtain $`\sigma _8=0.37\mathrm{\Omega }_0^{0.41+0.24\mathrm{\Omega }_0}`$ as the most probable value that results from the weak lensing analysis for the full SDSS/RASS sample. However, it turns out that changing the outer cluster mass density profile does not significantly mitigate the discrepancy between the values obtained for $`\sigma _8`$ using the two SDSS/RASS sub-samples, even for low values of $`\mathrm{\Omega }_0`$.

The discrepancy between the $`\sigma _8`$ results obtained for the two Sheldon et al. (2001) sub-samples could be alleviated if the mean projected mass of the high-luminosity sample is overestimated, or that of the low-luminosity sample underestimated. The first hypothesis is much more probable, and could be due to a contribution to the mean projected mass by filamentary material infalling into the clusters along the line of sight (Cen 1997; Metzler et al. 1999; Reblinsky & Bartelmann 1999; Metzler, White & Loken 2001). Alternatively, the mean X-ray luminosity of the high-luminosity sample may have been underestimated, or that of the low-luminosity sample overestimated. The latter could be caused by AGN contamination, but the initial study of Miller et al. (2002, in preparation) indicates that this is not a serious problem in the case of the SDSS/RASS cluster survey. Furthermore, the observed relationship between X-ray cluster luminosity and measured velocity dispersion (from the SDSS spectroscopic sample) for both the Sheldon et al. (2001) and the Miller et al. (2002, in prep.) samples are in excellent agreement with the $`L_\mathrm{X}`$$`\sigma _\mathrm{v}`$ relation of Mahdavi & Geller (2001).

Reiprich & Böhringer (2001) used ROSAT and ASCA X-ray data on 106 clusters to obtain the relation between X-ray luminosity in the \[0.1,2.4\] keV band and cluster mass in the form of $`M_{500}`$, with the cluster masses estimated assuming hydrostatic equilibrium (Finoguenov et al. 2001; Reiprich & Böhringer 2001). They found $`L_X=10^{18.59\pm 1.23}\times M_{500}^{1.575\pm 0.084}`$ where $`L_X`$ is in units of $`10^{40}\mathrm{erg}/\mathrm{s}`$ and $`M_{500}`$ in units of $`M_{}`$ (assuming $`h=0.5`$). Substituting the mean luminosities for the SDSS/RASS full sample and two sub-samples into this relation, one finds that the Sheldon et al. (2001) estimates for $`M_{500}`$ are well within the (extremely wide) range of possible values. Conversely, the relation between $`L_X`$ and $`M_{500}`$ can be estimated by combining the Sheldon et al. (2001) results with the shape of the cluster mass function, assumed well described by that of Jenkins et al. (2001), and the data on the cluster luminosity function from the REFLEX survey (Böhringer et al. 2001a,b). Assuming the $`L_XM_{500}`$ relation to be a power-law, we found that the normalization is essentially defined by the Sheldon et al. (2001) data, as expected, while the exponent is mainly governed by the relative shape of the mass and luminosity (cumulative) functions. Taking into account all the uncertainties involved in the normalization of the Jenkins et al. (2001) mass function by means of the Sheldon et al. (2001) data for the full SDSS/RASS sample, and those associated with the REFLEX luminosity function, we obtained through Monte Carlo simulations a 95 per cent confidence interval of \[1.6,3.7\] for the exponent of the $`L_X`$-$`M_{500}`$ relation. The allowed interval does not change significantly if the data for either of the two SDSS/RASS sub-samples is used instead to normalize the mass function. The low-luminosity sample yields \[1.8,3.2\], while from the high-luminosity sample we get \[1.5,3.1\]. In any case, surprisingly, the preferred value is close to 2.2, substantially higher than the $`1.58\pm 0.08`$ at $`1\sigma `$ found by Reiprich & Böhringer (2001). This analysis was performed for a flat universe with $`\mathrm{\Omega }_0=0.3`$, though very similar results were found for $`\mathrm{\Omega }_0=1`$. In order to confirm that the exponent of the $`L_XM_{500}`$ relation is only weakly determined by the normalization of the mass function, and thus by the Sheldon et al. (2001) data, we varied $`\sigma _8`$ between 0.5 and 1.2 (for $`\mathrm{\Omega }_0=0.3`$) and found that the preferred value for the exponent changes from 2.4 to 1.7. Therefore, assuming the Jenkins et al. (2001) mass function provides an accurate description of the cluster mass function, the discrepancy just found on the exponent of the $`L_XM_{500}`$ relation means that such a relation as obtained in Reiprich & Böhringer (2001) is at best only marginally consistent with the Sheldon et al. (2001) data taken together with the REFLEX luminosity function. And only if $`\sigma _8`$ is at the higher end of recent estimates (see e.g. Viana & Liddle 1999) can the Reiprich & Böhringer (2001) $`L_XM_{500}`$ relation be made consistent with the REFLEX luminosity function, within the context of the cosmological models discussed in this paper.

In calculating $`\sigma _8`$ we assumed that, on average, a galaxy cluster with the mean X-ray luminosity of the SDSS/RASS full sample has $`M_{500}`$ equal to the value estimated in Sheldon et al. (2001) for such a sample. This assumption is prone to several biases, one being that, assuming that there is some dispersion in the associated luminosity, the clusters of a given mass that preferentially end up in a sample selected in the manner of the SDSS/RASS are the most luminous ones. This bias leads to an underestimation of the correct mass corresponding to a given luminosity, and consequently to an underestimation of $`\sigma _8`$, though it is difficult to say by how much. The same effect takes place due to another, more subtle, type of bias, arising from the fact that the $`L_X`$$`M_{500}`$ relation is not linear. As a result, for a given sample the mean $`M_{500}`$ is not proportional to the mean $`L_X`$. Given that $`M_{500}L_X^n`$ with $`n<1`$, assuming that a cluster with the mean luminosity has a mass $`M_{500}`$ equal to the sample mean underestimates $`M_{500}`$ for such a cluster (because actually the mean $`M_{500}`$ is proportional to the mean of $`L_X^n`$). We have found that, for the luminosity dispersion of the clusters in the SDSS/RASS full sample, the most probable value of $`M_{500}`$ for a cluster with the mean $`L_X`$ may be underestimated by around 20 per cent. This percentage is robust to changes in the assumed value for the exponent $`n`$ (between 0.3 and 0.7) and to the possibility of dispersion in the $`L_X`$$`M_{500}`$ relation. This bias leads to a possible underestimation of $`\sigma _8`$ close to 10 per cent, almost independent of $`\mathrm{\Omega }_0`$. However, given that the mean $`M_{500}`$ for the SDSS/RASS sample is in fact not a mean of several independently calculated $`M_{500}`$, but the result of a mean shear profile, and that the mean $`L_X`$ is weighted by each cluster contribution to that profile, the above considerations may not be directly applicable to the case in hand. Note that because the dispersion in luminosities is smaller for either of the two SDSS/RASS sub-samples as compared to the full sample, the underestimation in $`\sigma _8`$ is smaller when it is estimated from the sub-samples, being closer to 5 per cent. Unfortunately, none of these biases seem to be able to significantly narrow the discrepancy between the $`\sigma _8`$ values derived using each SDSS/RASS sub-sample.

## 5 Conclusions

We have applied a new approach for constraining the normalization of the matter power spectrum, using the REFLEX X-ray luminosity function and the relation between cluster X-ray luminosity and mass obtained through weak lensing data from a preliminary small sample of SDSS/RASS clusters. We obtained $`\sigma _8`$ values significantly lower than other estimates based on cluster abundance data, with the exception of the recent results by Borgani et al. (2001), Reiprich & Böhringer (2001) and Seljak (2001). However, systematic biases may affect our analysis, given that barely consistent results are obtained when using subsets of the weak lensing data. This may be due to the small sample of clusters used, or an artifact of the techniques used in Sheldon et al. (2001) for co-adding clusters to produce an ensemble averaged weak lensing signal. In the process, we found that comparing the REFLEX luminosity function and the Jenkins et al. (2001) mass function implies that the relation between X-ray luminosity and cluster mass may be significantly steeper than previously thought.

The SDSS/RASS data set we have used will be dwarfed by the final SDSS/RASS catalogue (see Nichol et al. 2001), and surveys using the XMM-Newton satellite should supply much greater information on cluster luminosities (e.g. Romer et al. 2001). The prospect of considerably improving the constraint on $`\sigma _8`$ using this approach in the future is therefore great.

A.R.L. was supported in part by the Leverhulme Trust. We thank Hans Böhringer, Bob Mann and Kathy Romer for useful discussions and comments, Erin Sheldon for providing the mean X-ray luminosities of the SDSS/RASS samples in the ROSAT \[0.1,2.4\] keV band, and Vincent Eke for pointing out the possible bias in the $`\sigma _8`$ estimation due to the non-linear nature of the $`L_XM_{500}`$ relation.


text:
# Specific heat of MgB2 in a one- and a two-band model from first-principles calculations

## 1 Introduction

The nature of the superconducting state in MgB2 has been characterized by a broad range of experimental and theoretical methods and many basic properties have been unambiguously established since the discovery of the 40 K superconductor MgB2 by Nagamatsu and collaborators .

While electron-phonon coupling as the underlying pairing mechanism has been pinpointed by a large B isotope effect on $`T_c`$ proving B related vibrations to be essential  the nature of the order-parameter (viz. the superconducting gap) has remained a matter of debate. The order-parameter has been intensively investigated by tunnelling and point contact spectroscopy  as well as by high-resolution photoelectron spectroscopy . While these techniques show an energy gap in the quasiparticle spectrum most likely of $`s`$-wave symmetry the magnitude of the gap, $`\mathrm{\Delta }(0)`$, itself remained an open question: tunnelling experiments initially revealed a distribution of energy gaps with lower boundary 2$`\mathrm{\Delta }_1(0)`$/$`k_\mathrm{B}T_c`$ 1.1 and upper boundary 2$`\mathrm{\Delta }_2(0)`$/$`k_\mathrm{B}T_c`$ 4.5. These values are either considerably lower or distinctly larger than the weak coupling BCS value of 2$`\mathrm{\Delta }(0)`$/$`k_\mathrm{B}T_c=`$ 3.53 and these controversial findings have been discussed in terms of gap anisotropy or more recently attributed to the presence of two gaps or multiple gaps . The analysis of the electronic Raman continuum of MgB2 by Chen et al.  also pointed to the presence of two gaps with gap values within the limits indicated by the tunnelling experiments .

While these experiments employ surface sensitive techniques to determine the gap properties, evidence for multigap behaviour emerges also from methods like heat capacity or $`\mu SR`$ measurements probing true bulk properties . In the early heat capacity experiments the typical jump-like anomaly is seen at $`T_c`$ the magnitude $`\mathrm{\Delta }C_p`$/$`T_c`$ of which amounts at best to only about 70-80 % of the value 1.43 $`\gamma ^\mathrm{N}(T_c)`$ predicted by weak-coupling BCS theory . $`\gamma ^\mathrm{N}(T_c)`$ is the Sommerfeld constant in the normal state which was obtained from heat capacity measurements in high magnetic fields and which was determined to be 2.7-3 mJ mol-1 K-2. The shape of the heat capacity anomaly compares reasonably well with BCS-type behaviour assuming 2$`\mathrm{\Delta }(0)`$/$`k_\mathrm{B}T_c`$ = 3.53 with appropriately adjusted magnitude. An improved fit of the detailed temperature dependence of the heat capacity anomaly was obtained when calculating the heat capacity within the $`\alpha `$-model  assuming a BCS temperature dependence of the gap but with an increased ratio 2$`\mathrm{\Delta }(0)`$/$`k_\mathrm{B}T_c`$ = 4.2(2) . This result matches very well with the upper limit of the gap value consistently found in the tunnelling experiments and was suggested as an evidence that MgB2 is in the moderately strong coupling limit. More recently, the excess heat capacity observed close to $`T_c`$/4 by Bouquet et al.  and Wang et al.  has been attributed to a second smaller gap. Fits with a phenomenological two-gap model assuming that the heat capacity of MgB2 can be composed as a sum of the two individual heat capacities gave very good description with gap values of 2$`\mathrm{\Delta }_1(0)`$/$`k_\mathrm{B}T_c`$ = 1.2(1) and 2$`\mathrm{\Delta }_2(0)`$/$`k_\mathrm{B}T_c`$ 4 . Recent muon-spin-relaxation measurements of the magnetic penetration depth are consistent with a two-gap model .

Theoretically multigap superconductivity in MgB2 was first proposed by Shulga et al. to explain the behaviour of the upper critical magnetic field . Based on the electronic structure the existence of multiple gaps has been suggested by Liu et al. in order to explain the magnitude of $`T_c`$ . The electronic structure of MgB2 contains four Fermi surface sheets . Two of them with 2D character emerging from bonding $`\sigma `$ bands form small cylindrical Fermi surfaces around $`\mathrm{\Gamma }`$ \- A. The other two originating from bonding and antibonding $`\pi `$ bands have 3D character and form a tubular network. Liu et al. from first-principles calculations of the electron-phonon coupling conclude that the superconducting gap is different for the individual sheets and they obtain two different order parameters, a larger one on the 2D Fermi surface sheets and a second gap on the 3D Fermi surfaces, the latter was estimated to be approximately a factor of three reduced compared to the former .

In the present paper we calculate the specific heat capacity from the spectral Eliashberg function $`\alpha ^2(\omega )F(\omega )`$ first in the one-band model using the isotropic $`\alpha ^2(\omega )F(\omega )`$ as given by Kong et al.. Then we calculate the heat capacity in a two-band model by reducing the 16 Eliashberg functions $`\alpha _{ij}^2(\omega )F_{ij}(\omega )`$ appropriate for the four Fermi surface sheets into four Eliashberg functions corresponding to an effective-two-band model with a $`\sigma `$-band and $`\pi `$-band only. From the solution of the Eliashberg equations we obtain a superconducting gap ratio $`\mathrm{\Delta }_\sigma /\mathrm{\Delta }_\pi 2.63`$ in good agreement with the experimental data. The two-band model explains the reduced magnitude of the heat capacity anomaly at $`T_c`$ very well and also reproduces the experimental observed excess heat capacity at low temperatures.

## 2 Theory

### 2.1 One-Band Model

First we discuss the specific heat in the isotropic single band model with a strong (intermediate) electron-phonon interaction (EPI). In the normal state and in the adiabatic approximation the electronic contribution to the specific heat is determined from the Eliashberg function $`\alpha ^2(\omega )F(\omega )`$ by the expression

$`C_N^{el}(T)`$ $`=`$ $`(2/3)\pi ^2N(0)k_B^2T`$
$`\times \left[1+(6/\pi k_BT){\displaystyle _0^{\mathrm{}}}f(\omega /2\pi k_BT)\alpha ^2(\omega )F(\omega )\omega \right],`$
where $`N(0)`$ is a bare density of states per spin at the Fermi energy. The kernel $`f(x)`$ is expressed in terms of the derivatives of the digamma function $`\psi (x)`$

$$f(x)=x2x^2Im\psi ^{}(\mathrm{}x)x^3Re\psi ^{\prime \prime }(\mathrm{}x).$$
(2)
At low temperatures the specific heat has the well known asymptotic form: $`C_N^{el}(T0)=(1+\lambda )\gamma _0T`$, where $`\lambda =2_0^{\mathrm{}}𝑑\omega \omega ^1\alpha ^2(\omega )F(\omega )`$ is the electron-phonon coupling constant, and $`\gamma _0=2\pi ^2k_B^2N(0)/3`$ is the specific heat coefficient for noninteracting electrons. At higher temperatures the specific heat differs from this trivial expression (see, the discussion in reference ).

In the superconducting state an expression for the specific heat obtained by Bardeen and Stephen  which is based on an approximate sum rule has often been used. We shall however use an exact expression for the thermodynamical potential in the electron-phonon system which is based on the integration of the electronic Green’s function over the coupling constant

$$\mathrm{\Omega }=\mathrm{\Omega }_{\mathrm{el}}^{(0)}+\mathrm{\Omega }_{\mathrm{ph}}^{(0)}+T\underset{\omega _\mathrm{n}}{}\underset{0}{\overset{1}{}}\frac{x}{x}\mathrm{tr}\left[\widehat{\mathrm{\Sigma }}(x)\widehat{G}(x)\right]$$
(3)
where $`x`$ is dimensionless, $`\widehat{G}(x)`$ and $`\widehat{\mathrm{\Sigma }}(x)`$ are the exact electron Green’s function and the self-energy, respectively, for a coupling constant of $`x\lambda `$. The functions $`\mathrm{\Omega }_{\mathrm{el}}^{(0)}`$ and $`\mathrm{\Omega }_{\mathrm{ph}}^{(0)}`$ are the thermodynamic potentials for noninteracting electrons and noninteracting phonons, respectively. Some further arithmetics leads to the expression for the difference in free energies, $`F_\mathrm{N}`$ and $`F_\mathrm{S}`$, of the superconducting and normal state

$$\frac{F_NF_S}{\pi N(0)T}=\underset{n=\omega _c}{\overset{\omega _c}{}}\left\{\begin{array}{c}|\omega _n|(Z^N(\omega _n)1)\frac{2\omega _n^2[\left(Z^S(\omega _n)\right)^21]+\phi _n^2}{|\omega _n|+\sqrt{\omega _n^2\left(Z^S(\omega _n)\right)^2+\phi _n^2}}\\ +\frac{\omega _n^2Z^S(\omega _n)(Z^S(\omega _n)1)+\phi _n^2}{\sqrt{\omega _n^2\left(Z^S(\omega _n)\right)^2+\phi _n^2}}\end{array}\right\},$$
(4)
where $`Z(\omega _n)`$ is a normalization factor, $`\phi _n=\mathrm{\Delta }_n/Z(\omega _n)`$ is an order parameter, and $`\mathrm{\Delta }_n`$ is the gap function.

The specific heat at temperature, $`T`$, is calculated according to:

$$\mathrm{\Delta }C_{\mathrm{el}}(T)=T^2(F_NF_S)/T^2.$$
(5)
The specific heat jump $`\mathrm{\Delta }C_{\mathrm{el}}(T_c)`$ at $`T=T_c`$ is determined by the coefficient $`\beta =T_c\mathrm{\Delta }C_{\mathrm{el}}(T_c)/2`$ in $`F_NF_S=\beta t^2`$, where $`t=(T_cT)/T_c`$.

### 2.2 Two-Band Model

We have calculated the 16 Eliashberg functions $`\alpha _{ij}^2(\omega )F_{ij}(\omega )`$ where $`i`$ and $`j`$ label the four Fermi surface sheets and thereafter combined them into four corresponding to an effective two-band model which contains only a $`\sigma `$\- and $`\pi `$-band. Their respective densities of states at the Fermi energy have values of $`N_\sigma (0)=0.300`$ states/cell$``$eV and $`N_\pi (0)=0.410`$ states/cell$``$eV. Similar coupling constants $`\lambda _{\sigma \sigma },\lambda _{\pi \pi },\lambda _{\sigma \pi }`$, and $`\lambda _{\pi \sigma }`$ which are required for a two-band model were calculated earlier in reference. The procedure of reducing the 16 Eliashberg functions of the real 4 band system due to the 4 different Fermi surface sheets to an effective two-band model with only four coupling constants $`\lambda _{ij}`$ is an approximation which is based on the similarity of the two cylindrical and the two three-dimensional sheets of the Fermi surface requiring the same physical properties in both $`\sigma `$-bands or both $`\pi `$-bands. More details can be found elsewhere.

The four Eliashberg functions $`\alpha _{ij}^2(\omega )F_{ij}(\omega )`$ for the effective two-band model are shown in figure 1. The most significant contribution comes from the coupling of the bond stretching phonon modes to the $`\sigma `$-band. The coupling constants corresponding to the superconducting Eliashberg functions have been calculated to be: $`\lambda _{\sigma \sigma }=1.017,\lambda _{\pi \pi }=0.448,\lambda _{\sigma \pi }=0.213,`$ and $`\lambda _{\pi \sigma }=0.155`$. The small difference to the values given in reference  may be attributed to the different first-principles methods used in the calculation of the Eliashberg functions.

Besides the spectral functions we need to know the the Coulomb matrix element $`\mu _{ij}`$. With the help of the wavefunctions from our first-principles calculations we can approximately calculate the ratios for the $`\mu `$-matrix . The $`\sigma \sigma `$-, $`\pi \pi `$\- and $`\sigma \pi `$-values were in the ratio 2.23/2.48/1. This allows one to express $`\mu _{ij}^{}(\omega _c)`$ by these ratios and one single free parameter which is fixed to get the experimental $`T_c`$ of 39.4 K from the solution of the Eliashberg equations. The $`\mu ^{}(\omega _c)`$ matrix elements determined by this procedure are $`\mu _{\sigma \sigma }^{}(\omega _c)`$=0.210, $`\mu _{\sigma \pi }^{}(\omega _c)`$=0.095, $`\mu _{\pi \sigma }^{}(\omega _c)`$=0.069, and $`\mu _{\pi \pi }^{}(\omega _c)`$=0.172.

Using our calculated Eliashberg functions on the imaginary (Matsubara) axis together with the above matrix $`\mu _{ij}^{}(\omega _c)`$ we obtain the gap values $`\mathrm{\Delta }_\sigma =\underset{T0}{lim}\mathrm{\Delta }_\sigma (i\pi T)7.1`$meV, and $`\mathrm{\Delta }_\pi 2.7`$meV, which corresponds to $`2\mathrm{\Delta }_\sigma /k_BT_c`$=4.18 and $`2\mathrm{\Delta }_\pi /k_BT_c`$=1.59. The temperature dependence of the superconducting gaps is shown in figure 2. The filled circles (squares) display the gap for the 2D $`\sigma `$\- (3D $`\pi `$-) band. Due to interband coupling between the bands both gaps close at the same critical temperature. For a comparison also the BCS curve (line) is shown for a single gap (one-band model) which closes at $`T_c`$=39.4 K. The corresponding single BCS gap would be 6 meV.

The extension of equation 2.1 to the two-band model gives

$$C_{\mathrm{el}}^\mathrm{N}(T)=\frac{2\pi ^2}{3}N_{tot}(0)k_B^2T+4\pi k_B[N_\sigma (0)(I_{\sigma \sigma }+I_{\sigma \pi })+N_\pi (0)(I_{\pi \pi }+I_{\pi \sigma })]$$
(6)
where $`I_{ij}=_0^{\mathrm{}}f(\omega /2\pi k_BT)\alpha _{ij}^2(\omega )F_{ij}(\omega )𝑑\omega `$ ($`i,j=\pi ,\sigma `$), and the function $`f(x)`$ is given by equation 2.

The generalization of the superconducting free energy (4) to the two band model is straightforward and the heat capacity was obtained according to equation 5.

## 3 Comparison with Experiment

For the comparison with experiment we have selected data obtained by our group  and by Bouquet et al. The anomaly clearly visible at $`T_c`$ in the zero-field data is suppressed by a magnetic field of 9 Tesla in both experiments. In figure 3 we display the difference $`\mathrm{\Delta }C_p=C_p(0\mathrm{T}\mathrm{e}\mathrm{s}\mathrm{l}\mathrm{a})C_p(9\mathrm{T}\mathrm{e}\mathrm{s}\mathrm{l}\mathrm{a})`$. The anomalies at $`T_c`$ detected by both groups clearly have a different magnitude, the one described in reference  amounts to 133 mJ/mol K at $`T_c`$ and represents the largest specific heat capacity anomaly reported for MgB2 so far . The $`\mathrm{\Delta }C_p(T_c)`$ reported by our group is somewhat smaller, however, the shape of the anomalies close to $`T_c`$ is very similar for both samples. In fact, fitting the anomalies with the $`\alpha `$-model revealed an identical ratio $`2\mathrm{\Delta }/k_BT_c`$=4.2 with $`\mathrm{\Delta }`$=7 meV for both samples .

First we will try to discuss the experimental results in terms of a conventional one-band model. The specific heat in MgB2 was calculated using the isotropic spectral Eliashberg function $`\alpha ^2(\omega )F(\omega )`$ of Kong et al. . This function yields an electron-phonon coupling constant $`\lambda =0.87`$ and together with a Coulomb pseudopotential of $`\mu ^{}=0.1`$ yields $`T_c=40`$ K. The calculated specific heat at $`T_c`$ is $`\gamma ^\mathrm{N}(T_c)=1.94\gamma _0=3.24`$ mJ/mol K2 with $`\gamma _0`$=1.67 mJ/mol K2 from the band structure calculations of reference . The specific heat jump at $`T_c`$ equals $`\mathrm{\Delta }C`$ 196 mJ/mol K, which is a factor of 1.5-2 larger than the experimental values . It corresponds to $`\mathrm{\Delta }C/(\gamma ^\mathrm{N}(T_c)T_c)`$ 1.51 compared to the BCS value of 1.43. The difference $`\mathrm{\Delta }C_{\mathrm{el}}(T)=C_{\mathrm{el}}^S(T)C_{\mathrm{el}}^N(T)`$ is shown in figure 3 (dashed-dot line) in comparison with the experimental data. Not only the size of the jump disagrees with the experiment, but also the behaviour at low temperatures is different. The latter is connected with the fact that at low temperature equation (4) for a single band model leads to the standard exponential dependence $`C^ST^{3/2}\mathrm{exp}(\mathrm{\Delta }/T)`$, while the experimental data show a more complicated behaviour. Clearly there exists a discrepancy between experimental data and a theoretical one-band model.

The solid line in figure 3 represents the theoretical results for the two-band model as described above. The low temperature behaviour is in much better agreement with the experiment. The specific heat jump is now significantly reduced in comparison with a single band model and reproduces surprisingly well the experimental data of reference . With the data given above we obtain from our theoretical calculation an electronic heat capacity in the normal state of $`\gamma ^\mathrm{N}(T_c)=C_{\mathrm{el}}^\mathrm{N}(T_c)/T_c3.24`$ mJ/mol K2, the same value as for the one-band model. The absolute value of the specific heat jump in the two-band model is $`\mathrm{\Delta }C`$ 125 mJ/mol K, corresponding to $`\mathrm{\Delta }C/(\gamma ^\mathrm{N}(T_c)T_c)`$ 0.98 which is now smaller than the BCS value.

We would like to emphasize here that no fitting is involved in the theoretical calculations. The only free parameter which is in the Coulomb matrix elements is already determined by the experimental $`T_c`$ of 39.4 K.

One could expect that the difference between the theoretical results of the effective two-band model and our experimental data may be attributable to a different amount of impurities in the samples compared to the samples of reference . In the one-band model the critical temperature $`T_c`$ as well as the value and the temperature dependence of $`\mathrm{\Delta }C_p(T)`$ are not affected by non-magnetic impurities (Anderson theorem). This is in complete contrast to the situation for the two-band model, where both quantities are strongly dependent on interband impurity scattering. Interband impurity scattering leads to averaging of the gaps and thus to the increase of $`\mathrm{\Delta }C_p/\gamma ^\mathrm{N}(T_c)T_c`$ ratio. On the other hand, due to decrease of $`T_c`$, the specific heat jump only depends weakly on the scattering strength. In order to investigate the dependence of $`T_c`$ and of $`\mathrm{\Delta }C_p`$ on the interband impurity scattering we included the effect of interband impurities in the Eliashberg equations. The results show that even for rather strong impurity scattering $`1/2\tau =3\pi T_{c0}371`$ K, which leads to a drastic change of the critical temperature (decreasing to $`T_c`$=29.4 K) and strong averaging of the gaps, the specific heat jump remains practically unchanged $`\mathrm{\Delta }C_p`$ 120 mJ/mol K. This corresponds to a ratio $`\mathrm{\Delta }C_p/\gamma (T_c)T_c`$ 1.48, which is very close to the corresponding value of a single gap model. Therefore, interband impurity scattering can explain the change of $`T_c`$ in different samples, but is not responsible for the observed different values of the specific heat capacity anomaly at $`T_c`$.

We have shown that a complete theoretical calculation from first-principles using an effective two-band model can explain the major features in the specific heat measurement of MgB2 surprisingly well. The presented theoretical framework goes beyond a simple phenomenological two-gap model because interband effects are included explicitely and no fitting to experimental results has been performed. The reduced value of the heat capacity anomaly at $`T_c`$ as well as the low temperature behaviour are in excellent agreement with experimental results. The same first-principles approach using exactly the same Eliashberg functions and Coulomb matrix elements has been used in order to explain optical measurements  and tunnelling experiments  of the interesting superconductor MgB2.

JK would like to thank the Schloeßmann Foundation for financial support.

## References