text:
# Probability-Changing Cluster Algorithm for Two-Dimensional XY and Clock Models

## I Introduction

The two-dimensional (2D) XY model shows a unique phase transition, the Kosterlitz-Thouless (KT) transition . It does not have a true long-range order, but the correlation function decays as a power of the distance at all the temperatures below the KT transition point. José et al  studied the effect of the $`q`$-fold symmetry-breaking fields on the 2D XY model; this is essentially the same as treating the $`q`$-state clock model, where only the discrete values are allowed for the angle $`\theta _i`$ of the XY spins such that

$$\theta _i=2\pi p_i/q\mathrm{with}p_i=0,1,2,\mathrm{},q1.$$
(1)
In the limit $`q\mathrm{}`$ we get the XY model. It was shown that the 2D $`q`$-state clock model has two phase transitions of the KT type at $`T_1`$ and $`T_2`$ ($`T_1<T_2`$) for $`q>4`$. There is an intermediate XY-like phase between a low-temperature ordered phase ($`T<T_1`$) and a high-temperature disordered phase ($`T>T_2`$). It was predicted that the decay critical exponent $`\eta `$ varies from $`\eta =1/4`$ at $`T_2`$ to $`\eta =4/q^2`$ at $`T_1`$.

Most of the above theoretical analyses relied on the renormalization group argument, and they are not exact. There have been extensive numerical studies on the 2D classical XY model . In contrast, only a limited number of numerical works have been reported on the $`q`$-state clock model , and the accuracy was not good enough especially for the low-temperature phase transition. There have been no systematic studies to check the prediction of $`\eta (T_1)=4/q^2`$.

In numerical studies, efficient algorithms are important for getting the necessary information. The cluster update algorithms of the Monte Carlo simulation  are examples of such efforts, and they are useful for overcoming the problems of slow dynamics. Recently we proposed an effective cluster algorithm, which is called the probability-changing cluster (PCC) algorithm, of tuning the critical point automatically . It is an extension of the Swendsen-Wang algorithm , but we change the probability of cluster update (essentially, the temperature) depending on the observation whether clusters are percolating or not percolating. We showed the effectiveness of the PCC algorithm for the Potts models , determining the critical point and critical exponents for the second-order phase transition precisely with less numerical efforts. The PCC algorithm was also applied to the 2D site-diluted Ising model, where the crossover and the self-averaging properties were studied . The advantage of using the PCC algorithm in the study of random systems is as follows: The sample-dependent $`T_c(L)`$ for the finite system with the linear size $`L`$ is important for taking sample average, and the PCC algorithm is suitable for getting the sample-dependent $`T_c(L)`$.

There are a lot of interesting questions about the extension of the PCC algorithm. (i) Can the PCC algorithm be used for the problem of the vector order parameter, such as the XY model? (ii) Can it be applied to the analysis of the transition of the KT type? (iii) Can it work even if the system shows two or more phase transitions?

The purpose of the present paper is to answer these questions. We extend the PCC algorithm so as to treat systems with the vector order parameter. The rest of the paper is organized as follows. In Sec. II, we formulate the extension of the PCC algorithm for the vector order parameter. In Sec. III, the KT transition of the 2D XY model is studied with the finite-size scaling (FSS) analysis based on the KT form of the correlation length. In Sec. IV, we study the KT transitions of the 2D clock models. We investigate both phase transitions at $`T_1`$ and $`T_2`$ for the $`q=6,8,12`$ clock models. The summary and discussions are given in Sec. V.

## II PCC algorithm for vector order parameter

Our Hamiltonian is given by

$$=J\underset{i,j}{}𝑺_i𝑺_j,$$
(2)
where $`𝑺_i`$ is a planar unit vector, $`(\mathrm{cos}\theta _i,\mathrm{sin}\theta _i)`$, at site $`i`$; $`\theta _i`$ takes the value of $`[0,2\pi )`$ for the XY model, and the value given in Eq. (1) for the $`q`$-state clock model. The summation is taken over the nearest-neighbor pairs $`i,j`$.

In order to extend the PCC algorithm to systems with the vector order parameter, we use Wolff’s idea of the embedded cluster formalism . We project the vector $`𝑺_i`$ onto a randomly chosen unit vector $`𝒆_1`$ and another unit vector $`𝒆_2`$, perpendicular to $`𝒆_1`$, as

$$𝑺_i=𝒆_1\mathrm{cos}\varphi _i+𝒆_2\mathrm{sin}\varphi _i,$$
(3)
where $`\varphi _i`$ is the angle measured from the axis of the vector $`𝒆_1`$. Then, the Hamiltonian, Eq. (2), is rewritten as

$$=\underset{i,j}{}(J_{ij}^{(1)}ϵ_i^{(1)}ϵ_j^{(1)}+J_{ij}^{(2)}ϵ_i^{(2)}ϵ_j^{(2)})$$
(4)
with positive effective couplings

$$J_{ij}^{(1)}=J|\mathrm{cos}\varphi _i||\mathrm{cos}\varphi _j|,J_{ij}^{(2)}=J|\mathrm{sin}\varphi _i||\mathrm{sin}\varphi _j|$$
(5)
for two sets of Ising variables $`\{ϵ_i^{(1)}\}`$ and $`\{ϵ_i^{(2)}\}`$. Formally, we can restrict ourselves to the region $`[0,\pi /2)`$ for $`\{\varphi _i\}`$, and we write the partition function as

$`Z={\displaystyle _0^{\pi /2}}\{d\varphi _i\}`$ $`{\displaystyle \underset{\{ϵ_i^{(1)}=\pm 1\}}{}}\mathrm{exp}(\beta {\displaystyle \underset{i,j}{}}J_{ij}^{(1)}ϵ_i^{(1)}ϵ_j^{(1)})`$ (6)
$`\times `$ $`{\displaystyle \underset{\{ϵ_i^{(2)}=\pm 1\}}{}}\mathrm{exp}(\beta {\displaystyle \underset{i,j}{}}J_{ij}^{(2)}ϵ_i^{(2)}ϵ_j^{(2)})`$ (7)
with $`\beta =1/k_BT`$. Then, we can use the Kasteleyn-Fortuin (KF) cluster representation for the Ising spins . To make the KF cluster, we connect the bonds of parallel Ising spins with the probability

$$p_{ij}^{(1,2)}=1\mathrm{exp}(2\beta J_{ij}^{(1,2)}).$$
(8)
In the PCC algorithm , the cluster representation of the Ising model is used in two ways. First, we flip all the spins on any KF cluster to one of two states, that is, $`+1`$ or $`1`$, as in the Swendsen-Wang algorithm . Second, we change the KF probability, Eq. (8), depending on the observation whether clusters are percolating or not. It is based on the fact that the spin-spin correlation function $`G(𝒓_i𝒓_j)`$ becomes nonzero for $`|𝒓_i𝒓_j|\mathrm{}`$ at the same point as the percolation threshold. For the XY model in the embedded cluster formalism, the spin-spin correlation function is written as

$`G(𝒓_i𝒓_j)`$ $`=`$ $`𝑺_i𝑺_j=|\mathrm{cos}\varphi _i||\mathrm{cos}\varphi _j|ϵ_i^{(1)}ϵ_j^{(1)}`$ (9)
$`+|\mathrm{sin}\varphi _i||\mathrm{sin}\varphi _j|ϵ_i^{(2)}ϵ_j^{(2)}`$ (10)
$`=`$ $`A_{ij}\mathrm{\Theta }^{(1)}(𝒓_i,𝒓_j)+B_{ij}\mathrm{\Theta }^{(2)}(𝒓_i,𝒓_j),`$ (11)
where $`\mathrm{}`$ represent the thermal average. The function $`\mathrm{\Theta }(𝒓_i,𝒓_j)`$ is equal to 1 (0) if the sites $`i`$ and $`j`$ belong to the same (different) cluster, and $`A_{ij}`$ and $`B_{ij}`$ are some constants. Thus, the system is regarded as percolating if $`ϵ^{(1)}`$ or $`ϵ^{(2)}`$ Ising spins are percolating. When treating the cluster spin update, one may consider Ising spins of a single type projected onto a randomly chosen axis, as in the original Wolff’s proposal . However, we should consider Ising spins of both types perpendicular to each other for checking the percolation.

The procedure of Monte Carlo spin update is as follows: (i) Start from some spin configuration and some value of $`\beta `$. (ii) Choose a unit vector $`𝒆_1`$ randomly. (iii) Construct the KF clusters for $`ϵ^{(1)}`$ and $`ϵ^{(2)}`$ using the probability, Eq. (8), and check whether the system is percolating or not. Flip all the spins on any KF cluster to $`+1`$ or $`1`$ for both $`ϵ^{(1)}`$ and $`ϵ^{(2)}`$ Ising spins. (iv) If the system is percolating (not percolating), decrease (increase) $`\beta `$ by $`\mathrm{\Delta }\beta (>0)`$. (v) Go back to the process ii.

As $`\mathrm{\Delta }\beta `$ becomes small, the distribution of $`\beta `$ becomes a sharp Gaussian distribution around the mean value $`\beta _c(L)`$, which depends on the system size $`L`$. We approach the canonical ensemble in this limit, and the existence probability $`E_p`$, the probability that the system percolates, becomes 1/2 at $`\beta _c(L)`$.

## III XY model

We have made simulations for the classical XY model on the square lattice with the system sizes $`L`$ =8, 16, 32, 64, 128, 256, and 512. After 20,000 Monte Carlo sweeps of determining $`\beta _c(L)`$ with gradually reducing $`\mathrm{\Delta }\beta `$, we have made 10,000 Monte Carlo sweeps to take the thermal average; we have made 100 runs for each size to get better statistics and to evaluate the statistical errors. As for the criterion to determine percolating, we have employed the topological rule  in the present study. The topological rule is that some cluster winds around the system in at least one of the $`D`$ directions in $`D`$-dimensional systems.

Let us start with the size dependence of the transition temperature. We use the FSS analysis based on the KT form of the correlation length,

$$\xi \mathrm{exp}(c/\sqrt{t})$$
(12)
with $`t=(TT_{\mathrm{KT}})/T_{\mathrm{KT}}`$. Using the PCC algorithm, we locate the temperature $`T_{\mathrm{KT}}(L)=1/k_B\beta _c(L)`$ that the existence probability $`E_p`$ is 1/2. Then, using the FSS form of $`E_p`$, that is, $`E_p=E_p(\xi /L)`$, we have the relation

$$T_{\mathrm{KT}}(L)=T_{\mathrm{KT}}+\frac{c^2T_{\mathrm{KT}}}{(\mathrm{ln}bL)^2}.$$
(13)
We plot $`T_{\mathrm{KT}}(L)`$ as a function of $`l^2`$ with $`l=\mathrm{ln}bL`$ for the best fitted parameters in Fig. 1. We represent the temperature in units of $`J/k_B`$. The error bars are smaller than the size of marks. Our estimate of $`T_{\mathrm{KT}}`$ is 0.8933(6); the number in the parentheses denotes the uncertainty in the last digits. We have estimated the uncertainty by the $`\chi ^2`$ test of the data for 100 samples. This value is consistent with the estimates of recent studies; 0.89213(10) by the Monte Carlo simulation , and 0.894 by the short-time dynamics . The constant $`c`$, in Eq. (13), is estimated as $`c`$=1.73(2).

Let us consider the magnetization $`m^2(L)`$ at $`T_{\mathrm{KT}}(L)`$ to discuss the critical exponent $`\eta `$. In Fig. 2, we plot $`m^2(L)`$ as a function of $`L`$ in logarithmic scale. We expect the FSS of the form $`m^2(L)L^\eta `$, but there are small corrections. The importance of the multiplicative logarithmic corrections were pointed out . Using the form

$$m^2(L)=AL^\eta (\mathrm{ln}b^{}L)^{2r},$$
(14)
we obtain $`\eta =0.243(5)`$ and $`r=0.038(5)`$. We show the fitting curve obtained by using Eq. (14) in Fig. 2. This value of $`\eta `$ is a little bit smaller than the theoretical prediction, 1/4 (=0.25). Our logarithmic-correction exponent $`r`$ is compatible with Janke’s result $`r=0.0560(17)`$ for thermodynamic data , but different from the theoretical prediction $`r=1/16`$ .

## IV Clock model

Next turn to the $`q`$-state clock model. Because of the reflection symmetry, we confine ourselves to the case of even $`q`$. Then, the same procedure can be used as the XY model. One thing we should have in mind is that the axis of the vector $`𝒆_1`$ should be chosen from one of $`q`$ directions in Eq. (1) or the middle of two of them. We plot the high-temperature transition temperature $`T_2(L)`$ of the 6-state clock model as a function of $`l^2`$ in Fig. 3. The estimate of $`T_2`$ is 0.9008(6), which is more precise than the previous estimates; 0.92(1)  and 0.90 . The plot of $`m^2(L)`$ at $`T_2(L)`$ as a function of $`L`$ for the 6-state clock model is given in Fig. 4. The estimate of $`\eta `$ is 0.243(5) by the analysis of the multiplicative logarithmic corrections, Eq. (14), and the exponent $`r`$ is estimated as 0.037(5).

For the second-order transition, the curves of the existence probability $`E_p`$ of different sizes cross at $`T_c`$ as far as the corrections to FSS are negligible; this is the same as the behavior of the Binder ratio . For the KT transition, however, $`T_{\mathrm{KT}}`$ is not the crossing point but the spray-out point. Therefore, $`T_2`$ can be searched only from the high-temperature side, and $`T_1`$ only from the low-temperature side. The value of $`E_p`$ at $`T_1`$ is close to 1. In principle, we can use the same procedure as the study of $`T_2`$; we may change the setting value of $`E_p`$, 1/2, to a higher one by introducing a biased random walk. However, it is difficult to resolve the size dependence for lower temperatures. Therefore, we employ a slightly different approach for the analysis of the phase transition at $`T_1`$. When judging whether clusters are percolating or not, we consider another type of clusters. Instead of choosing the vector $`𝒆_1`$ randomly in Eq. (3), we choose the vector $`𝒆_1`$ as

$$𝑴=|𝑴|𝒆_1$$
(15)
with $`𝑴=_i𝑺_i`$, or more generally we may choose $`𝒆_1`$ such that

$$𝑴=|𝑴|(𝒆_1\mathrm{cos}\varphi +𝒆_2\mathrm{sin}\varphi )$$
(16)
with some fixed angle $`\varphi `$. With this choice, the existence probability for the percolation of only $`ϵ^{(1)}`$ (or $`ϵ^{(2)}`$) Ising spins holds the same FSS property as the total $`E_p`$. As a result, we can control the value of $`E_p^{(1)}`$ at $`T_1`$ so as to apply the FSS analysis easily with an appropriate $`\varphi `$.

The low-temperature transition temperature $`T_1(L)`$ of the 6-state clock model obtained by the above modified approach is plotted as a function of $`l^2`$ also in Fig. 3. As the angle $`\varphi `$ in Eq. (16), we have used $`\pi /3`$. Our estimate of $`T_1`$ is 0.7014(11), which is more precise again than the previous estimates; 0.68(2)  and 0.75 . The plot of $`m^2(L)`$ at $`T_1(L)`$ for the 6-state clock model is also given in Fig. 4. The estimate of $`\eta `$ is 0.113(3) by the analysis of the multiplicative logarithmic corrections, Eq. (14), and the exponent $`r`$ is estimated as 0.017(4). The previous estimates of $`\eta `$ are 0.100  and 0.15 .

We have also made simulations for $`q=8`$ and $`q=12`$. The estimates of the transition temperatures $`T_1,T_2`$ and those of $`\eta (T_1)`$ and $`\eta (T_2)`$ for $`q`$=6, 8, 12, and $`\mathrm{}`$ (the XY model) are tabulated in Table I. The $`1/q^2`$-dependence of transition temperatures and exponents are shown in Fig. 5. There, the exact results for $`q=4`$ are also given; that is, the Ising singularity at $`T_c=[\mathrm{ln}(\sqrt{2}+1)]^1=1.1346`$ with $`\eta =1/4`$. The transition temperature $`T_1`$ becomes smoothly lower with larger $`q`$; in the lowest order we find that $`T_11/q^2`$, which is consistent with the theoretical prediction . The critical exponent $`\eta `$ at $`T_2`$ is a universal constant, and compatible with the theoretical prediction $`\eta =1/4`$. The estimates of the critical exponent $`\eta `$ at $`T_1`$ remarkably coincide with the theoretical prediction $`\eta =4/q^2`$; 1/9=0.111 for $`q`$=6, 1/16=0.0625 for $`q`$=8, and 1/36=0.0278 for $`q`$=12. This is the first systematic report of confirming the theoretical prediction.

## V Summary and discussions

To summarize, we have extended the PCC algorithm  to the study of the XY and clock models. Wolff’s idea of the embedded cluster formalism  is used for treating the system with the vector order parameter. The KT transitions of the 2D XY and clock models are studied by using the FSS analysis based on the KT form of the correlation length. For dealing with the low-temperature transition temperature, $`T_1`$, we have employed a slightly modified algorithm. Investigating the $`q=6,8,12`$ clock models, we have systematically confirmed the prediction of $`\eta (T_1)=4/q^2`$. We have shown that small logarithmic corrections are present in the KT transitions. The sign of the logarithmic-correction exponent $`r`$ is positive for all cases of the XY model and the clock models at both $`T_1`$ and $`T_2`$, which is compatible with Janke’s result , but different from the theoretical prediction $`r=1/16`$ . The present precise numerical results may stimulate the refined renormalization-group study of the KT transitions.

In the previous numerical studies of the KT transitions, one might resort to a big scale calculation using an extensive computer resource , or one might use some special boundary conditions . It is due to the subtlety of the KT phase transitions; that is, $`T_{\mathrm{KT}}`$ is not the crossing point but the spray-out point of the existence probability or the Binder parameter. Moreover, the low-temperature transition $`T_1`$ for the clock model is difficult to study because the system is nearly ordered; it is more difficult for larger $`q`$. We should stress that using the present efficient method of numerical simulation, we have successfully made a systematic study of the XY and clock models with much less efforts.

Our formalism of the vector order parameter is easily extended to the general $`O(n)`$ model, where the percolation of $`n`$ types of Ising spins will be considered. Then, more problems of interest can be studied by the PCC algorithm. The PCC algorithm can be also applied to the quantum Monte Carlo simulation with the cluster algorithm . It will be interesting to compare the present result with the quantum XY model.

## Acknowledgments

We thank N. Kawashima, H. Otsuka, M. Itakura, and Y. Ozeki for valuable discussions. Thanks are also due to M. Creutz for bringing our attention to the clock model having two phase transitions. The computation in this work has been done using the facilities of the Supercomputer Center, Institute for Solid State Physics, University of Tokyo. This work was supported by a Grant-in-Aid for Scientific Research from the Japan Society for the Promotion of Science.


format: text

text:
# Thermoelectric power of MgB2-xBex

## I Introduction

Since the recent discovery of superconductivity in MgB2, both theoretical and experimental efforts have been made to understand its structural and electrical properties. Most theoretical works suggest that coupling of the $`\sigma `$-hole with B-plane phonon is the key ingredient of the superconductivity. Boron isotope effect and other experimental data showed that the material is in the intermediate or strong BCS coupling regiem ($`\lambda =0.71`$).

Chemical substitution experiment has drawn much attention due to the possibility of obtaining higher superconducting transition temperature $`T_c.`$ In general, by replacing Mg- or B-sites with other elements, physical parameters such as lattice constants and carrier density change. Therefore, study of $`T_c`$ shift in well-controlled substitution samples provides a chance to understand the superconductivity in detail. For the Mg-site substitution, several compounds such as Mg1-xAlxB2 have been studied. For the B-site, C- and Be-substitutions have been reported.

It is important to know how these substitutions change the carrier density and how the change is related with $`T_c`$ shift. Thermoelectric power (TEP) measurement is a useful probe of the carrier density. For example, TEP experiment on Mg1-xAlxB2 showed that the Al-substitution dopes electrons. Theoretical calculation on Mg1-xAlxB2 suggested that the doping occurs largely in the $`\sigma `$-hole pocket. In the present paper, we report our results of TEP measurement and theoretical band calculation on MgB2-xBex samples where B is substituted with Be. In an earlier paper, we showed that in this compound, the MgB2 phase is maintained up to $`x=0.6.`$ Also the lattice constants and $`T_c`$ changed systematically with Be-substitution. We find that the TEP decreases with $`x`$, opposite to that in Mg1-xAlxB2, which suggests that hole is doped. Also, changes in the band structure indicate that the hole doping occurs mostly in the $`\sigma `$-band. Thus, Be-substitution brings about hole-doping into the boron plane. We consider the consequence of this effect on the transition temperature and arue that the hole density change plays only a second role in the $`T_c`$ shift.

Meanwhile, it is well known now that the temperature dependence of TEP is one of the unconventional features of MgB2: At low $`T,`$ TEP is linear in $`T`$ which is normal for most metals, but crosses over to a sublinear behavior as $`T`$ increases. There have been many interpretations of this phenomenon.  We find that the crossover exhibits a systematic change with the Be-substitution. Furthur, we show that TEP in the whole $`T`$ region is explained by a single model function in which the electron-phonon interaction is explicitely taken into account.

## II Experimental

Polycrystalline samples used in this experiment were synthesized by a powder metallurgical technique using a high pressure furnace. Starting materials are fine powders (-325 mesh) of Mg (99.8%, Alfa Aesar), amorphous B (99.99%, Alfa Aesar), and Be (99.+%, Alfa Aesar). Stoichiometric amounts of powders are mixed and pelletized. The pellets are placed in a tungsten vessel with a close-fitting cap, reacted two hour at 850 C under 20 atm. of high purity argon atmosphere.

In MgB2-xBe$`_x,`$ the MgB2 phase is maintained up to $`x=0.6`$. Their structural and superconducting properties were reported elsewhere. We found that as a function of Be-substitution, lattice parameters show decreasing $`a`$\- and increasing $`c`$-values, and transition temperature $`T_c`$ decreases as summarized in Table I.

For TEP measurements, bar-shaped samples (with dimensions of $``$ 4$`\times `$0.5$`\times `$0.1 $`mm^3`$) were mounted on two resistive heaters. Gold wires were used for thermoelectric potential leads. Chromel-constantan thermocouples were used for the temperature gradient measurement. Sample ends and thermocouple beads were glued to the heater blocks by Stycast epoxy. In our measurement, thermopower from the contact wires was carefully calibrated.

## III Result and Discussion

Figure 1 shows thermoelectric power $`S(T)`$ of MgB2-xBex samples. For MgB2 ($`x=0`$), $`S(T)`$ is linear in $`T`$ at low temperature above $`T_c`$. At higher $`T`$ ($`>150`$ K), it crosses over to a sub-linear behavior. These features together with the large jump at $`T_c`$ ($``$ 1.4 $`\mu `$V/K) are consistent with the earlier reports. As boron is substituted with Be, the linear slope decreases. Also the crossover temperature is lowered. At $`x=0.6`$, $`S(T)`$ changes sign at a low temperature.

Thermal diffusion of carriers gives rise to the linear-in-$`T`$ behavior in $`S(T)`$ at low temperatures. Kinetic transport theory shows that $`S(T)=(\pi ^2k_B^2T/3e)\sigma ^{^{}}(E)/\sigma (E)`$, where the dc-conductivity $`\sigma (E)`$ and its derivative $`\sigma ^{^{}}(E)`$ are calculated at the Fermi energy.  In a single parabolic band system, it is approximated to the Mott formula:

$$S(T)=X_bT=\pm \left(\frac{\pi ^2k_B^2}{3eE_F}\right)T,$$
(1)
where $`k_B`$ is the Boltzman constant, $`\pm `$corresponds to the carrier sign, and $`E_F`$ is the Fermi energy relative to the band maximum (or minimum). However, MgB2 is a multi band system where $`\sigma `$\- and $`\pi `$-bands coexist. Also the Fermi surface is not spherical. Thus, Eq. (1) can be applied only approximately here. The positive value of $`X_b`$ in pristine MgB2 suggests that the dominant carrier is hole. The decrease of $`X_b`$ with $`x`$ suggests that the hole density increases. To obtain $`X_b`$, we fit the linear part of $`S(T)`$ with Eq. (1) as shown by the dashed lines. Since the data are not extrapolated to zero, vertical shifts were needed in the fit.

To understand the behavior of $`X_b`$ more quantitatively, we calculated band structures of MgB2-xBex as shown in Fig. 2. Here we used the local density approximation (LDA) with the linearized augmented plane wave method. To account for the random substitution of Be, the virtual crystal approximation (VCA) was employed. Details of the calculation method were described by Mehl et al. In current work, the experimental lattice parameters in Table I were used. Four bands contribute to the Fermi surfaces: two $`\sigma `$-bands with B $`p_{x,y}`$ character give the 2D hole-type cylinders around the $`\mathrm{\Gamma }`$-$`A`$ line and two $`\pi `$-bands form the 3D honeycomb tubular networks. The latter consists of one electron-type at the $`H`$-point and another hole-type at the $`K`$-point. Note that the most prominent change with $`x`$ is the hole increase in the $`\sigma `$-bands. Change in the $`\pi `$-bands is relatively small. These results tell us that the Be-substitution dopes hole into the $`\sigma `$-bands.

It is interesting to note that $`T_c`$ decreases with the hole-doping, similarly to the electron-doped Mg1-xAlxB2. Theoretical studies show that the superconductivity arise from the $`\sigma `$-band holes coupled with the B-plane phonons. $`T_c=39`$ K in MgB2 can be produced from the McMillan formular using the coupling constant $`\lambda =1.01`$ and the Coulomb pseudopotential $`\mu ^{}=0.13`$. Here $`\lambda `$ is proportional to the density of state (DOS) of the $`\sigma `$-bands at the Fermi energy, $`N_h(E_F).`$ From our band calculation, we find that $`N_h(E_F)`$ =0.22 (eV-1 per cell) at $`x=0`$ and $`N_h(E_F)`$ =0.24 (eV-1 per cell) at $`x=0.6.`$ (the $`\sigma `$-bands are highly 2D-like and the DOS increases only slightly with hole doping.) Then $`\lambda `$ will increase proportionally to become 1.09 at $`x=0.6,`$ if we assume the other parameters determining $`\lambda `$ do not change. (see for example Eq. (3) of Ref. 3) This yields $`T_c=45`$ K, which is in sharp contrast with the observed $`T_c`$ decrease. This suggests that the other parameters change significantly with the substitution and their effects overcome the $`N_h(E_F)`$ contribution. In another paper, we dealt with this issue and showed that the lattice constant change is the primary cause of the reduced superconductivity.

Now we consider the unusual behavior of $`S(T)`$, i.e., the deviation from the linear dependence at high $`T`$. In previous works, it has been attributed to the minor carrier contribution, the thermally activated transport, and to the phonon-drag effect. Here, we consider effect of electron-phonon interaction on $`S(T).`$ According to Kaiser, the interaction contributes to enhance the TEP through modifying the carrier mass and thus the thermal diffusion. Taking this effect into account, Eq. (1) is rewritten as

$$S(T)=[1+\lambda \overline{\lambda }_s(T)]X_bT,$$
(2)
where $`X_b`$ is the slope in Eq. (1), $`\lambda `$ is the electron-phonon coupling constant, and $`\overline{\lambda }_s(T)`$ is a function which represents the $`T`$-dependent thermopower enhancement:

$$\overline{\lambda }_s(T)=_0^{\mathrm{}}𝑑\omega \frac{\alpha ^2F(\omega )}{\omega }G_s(\frac{\mathrm{}\omega }{k_BT}).$$
(3)
Here, the normalized Eliashberg function $`\alpha ^2F(\omega )`$ consists of the phonon density of states $`F(\omega )`$ and the coupling constant $`\alpha `$. $`G_s(\mathrm{}\omega /k_BT)`$ is a function associated with thermal population of phonons. For MgB$`_2,`$we calculated $`\overline{\lambda }_s(T)`$ using $`\alpha ^2F(\omega )`$ reported by Liu et al. and fit the data with Eq. (2). Here we used $`\lambda `$ and $`X_b`$ as fitting parameters. For $`x>0,`$ the Be-substitution into the B-plane will change, probably significantly, the phonon structure. Thus $`\alpha `$ and $`F(\omega )`$ will depend on $`x`$. Since they are not known, we took the values of MgB2 in calculating $`\overline{\lambda }_s(T)`$. Thus $`\lambda `$ and $`X_b`$ we estimate for $`x>0`$ samples are under large uncertainties.

Figure 3(a) shows the fit for MgB2 (solid line). The bare diffusion part ($`X_bT`$) and the enhancement part ($`\lambda \overline{\lambda }_s(T)\times X_bT`$) are represented with dashed and dash-dotted line, respectively. The latter has a broad maximum at $`T215`$ K. Inset shows calculated behavior at higher temperature. Fig. 3(b) shows the fitting results for MgB2-xBex. Note that the fit is reasonable except the small deviation for $`x=0.3`$.

In Fig. 4, we summarize the linear slope $`X_b`$ obtained from our analyses. The fitting results using Eq. (1) and Eq. (2) are shown with the filled-circles and the triangles, respectively. Note that $`X_b`$ from the modified model is smaller than that from the bare diffusion model. This is due to the enhancement effect contained in Eq. (2). Also, we estimated $`X_b`$ from the band calculation (the dash-dotted line). Here the $`E_F`$ in Eq. (1) was taken from the $`\sigma `$-hole bands, assuming contributions from the other bands are negligible. The band calculation result is closer to the modified model result at low doping region, $`x0.3`$, which supports the importance of the electron-phonon interaction effect. At $`x=0.3`$, $`X_b`$ from the two fits exhibit a sudden drop. This drop may be related to the observed structural change in the same composition. The incomplete agreement between the fit result and the band calculation result may come from complex effects not included in this work such as the multi-band contributions and anisotropic transport.

Regarding the electron-phonon coupling constant (EPC), we obtain $`\lambda =0.90`$ for MgB2. This is in good agreement with the earlier reports of 0.7 $``$ 1.0. For $`x>0`$, $`\lambda `$ increases to 0.98($`x=0.2`$), 1.31(0.3), 1.34(0.4), 1.47(0.6). This result is quite unusual because, as $`T_c`$ decreases with $`x,`$ $`\lambda `$ is expected to decrease. Recently, evidences show that MgB2 has two gaps. In this case, EPC from transport measurement (= $`\lambda _{tr})`$ is different from the EPC which determines the superconducting $`T_c`$ (= $`\lambda _{sc}`$). Thus, the increase of TEP $`\lambda `$ (= $`\lambda _{tr}`$) does not necessarily contradict with the $`T_c`$ decrease. One should also keep in mind that the increase of TEP $`\lambda `$ may be simply an errorneous effect that arise from the uncertainties in $`\alpha `$ and $`F(\omega )`$ for $`x>0`$ we mentioned above.

Now, let us consider the sign change in $`S(T)`$ of $`x=0.6`$. Sign change in TEP is widely observed in many alloyed systems, for example, Ag-Au alloy,  YBa2Cu3O7-δ, NbNx, , etc. In YBa2Cu3O7-δ, the change is observed as oxygen deficiency $`\delta `$ increases. In NbNx, $`S(T)`$ is composed of the diffusive (positive in sign) and phonon-drag terms (negative in sign). At high $`T`$, the former is dominant while the latter prevails at low $`T`$. In the intermediate $`T`$ , sign change occurs. It is tempting to interpret our observation similarly: the negative $`S(T)`$ may correspond to the phonon-drag effect. However, note that in the pristine MgB2, the phonon-drag feature is not observed. Furthur, the feature, if any, should be suppressed with Be-substitution because the phonon-drag generally disappears as randomness is increased. Origin of the sign change is thus remains for future study.

## IV Conclusion

¿From the TEP measurement and band structure calculation on MgB2-xBex ($`x=0`$, 0.2, 0.3, 0.4, and 0.6), we found that the hole density increases with $`x`$ in the $`\sigma `$-bands . Thus, the Be-substitution dopes hole into the boron plane. The fact that $`T_c`$ shifts in the same direction (lowering) as the electron-doped case suggests that carrier doping is not the primary route to control the transition temperature in MgB$`_2.`$ This result is consistent with the 2D nature of the $`\sigma `$-bands. Furthur, we showed that the anomalous behavior of TEP at high temperature can be explained.by taking the electron-phonon interaction effect into account.

## Acknowledgments

This work was supported by KRF-99-041-D00185 and by the KOSEF through the CSCMR.


text:
# Non-rigidity degrees of root lattices and their duals

## 1 Introduction

Voronoi \[Vo1908\] defined the partition of the cone $`𝒫_n`$ of positive semidefinite $`n`$-ary quadratic forms into $`L`$-type domains, that we call here $`L`$-domains. Forms of the same $`L`$-domain correspond to lattices that determine affinely equivalent Voronoi partitions of $`𝐑^n`$, i.e. partition into Voronoi polytopes. The partition of $`𝐑^n`$, which is (both, combinatorially and affinely) dual to a Voronoi partition, is called Delaunay partition and consists of Delaunay polytopes. In other words, two lattices have the same $`L`$-type if and only if the face posets of their Voronoi polytopes are isomorphic. (The face poset of a polytope $`P`$ is the set of all faces of $`P`$ of all dimensions ordered by inclusion.)

If the $`L`$-type of a lattice changes, then either some Delaunay polytopes are glued into a new Delaunay polytope, or a Delaunay polytope is partitioned into several new Delaunay polytopes. Recall that the center of a Delaunay polytope is a vertex of a Voronoi polytope. Hence if the $`L`$-type of a lattice changes, then for each Voronoi polytope either some vertices are glued into one vertex, or a vertex splits into several new vertices.

Voronoi proved that each $`L`$-domain is an open polyhedral cone of the dimension $`k`$, $`1kN`$, where $`N=\frac{n(n+1)}{2}`$ is the dimension of $`𝒫_n`$, i.e. dimension of the space of coefficients of $`f𝒫_n`$. An $`L`$-type having an $`N`$-dimensional domain, is called generic. Otherwise, $`L`$-type is called special.

In \[BG01\], a notion of non-rigidity degree of a form $`f`$ and the corresponding lattice $`L(f)`$ is introduced. It is denoted by nrd$`f`$ and is equal to dimension of the $`L`$-domain containing $`f`$. It is shown in \[BG01\] that nrd$`f`$ is equal to corank of a system of equalities connecting the norms of minimal vectors of cosets $`2L`$ in $`L`$. In fact, nrd is the number of degrees of freedom, such that the Delaunay partition (more precisely, its star) has, when one deforms it affinely, so that the result remains an Delaunay partition. Clearly, the maximal nrd is $`\left(\genfrac{}{}{0pt}{}{n+1}{2}\right)`$ and it is realized only by simplicial Delaunay partitions.

So, $`1\mathrm{nrd}fN`$, and nrd$`f=N`$ if $`f`$ belongs to a generic $`L`$-domain. A form $`f`$ and the corresponding lattice $`L(f)`$ are called rigid if nrd$`f=1`$. This name was used because any, distinct from a scaling, affine transformation of a rigid lattice changes its $`L`$-type. (Sometimes, a rigid form is called an edge form, since it lies on an extreme ray of the closure of an $`L`$-domain.) Clearly, any 1-dimensional lattice is trivially rigid. In \[BG01\], were given 7 examples of rigid lattices of dimension 5 and shown that $`D_4`$ is unique such a lattice of dimension $`n`$, $`2n4`$.

The set of all $`L`$-domains is partitioned into classes of unimodularly equivalent domains. For $`n3`$, there is only one class, i.e. only one generic $`L`$-type. For $`n=4`$, there are 3 generic $`L`$-types.

The $`L`$-domain $`𝒜_n`$ of the lattice $`A_n^{}`$ which is dual to the root lattice $`A_n`$ is well known. $`𝒜_n`$ has dimension $`N`$, i.e. it is generic. It is the unique $`L`$-domain such that all the extreme rays of its closure span forms of rank 1. All these facts were known to Voronoi. He called the domain $`𝒜_n`$ the first type domain and one of the forms of $`A_n^{}`$ the principal form of the first type.

The $`L`$-domain of the lattice $`A_n`$, $`n2`$, is a simplicial $`(n+1)`$-dimensional cone; it is described in \[BG01\]. So, nrd$`A_n=n+1`$ for $`n2`$. Also nrd$`Z_n=n`$ for $`n1`$ and $`A_1,D_2,D_3`$ are scalings of $`Z_1,Z_2,A_3`$, respectively. So, nrd$`(A_1=A_1^{})=1`$, nrd$`(D_2=D_2^{})=2`$ and nrd$`D_3=4`$, nrd$`D_3^{}=6`$.

It is proved in \[BG01\] that the lattice $`D_n`$ is rigid for $`n4`$.

The lattices $`E_6`$, $`E_6^{}`$, $`E_7`$, $`E_7^{}`$ and $`E_8=E_8^{}`$ are rigid, i.e. their $`L`$-domains are one-dimensional. The rigidity of root lattices $`E_6`$, $`E_7`$ and $`E_8`$ is shown in \[BG01\]. The rigidity of the lattice $`E_6^{}`$ was proved independently by Engel and Erdahl (personal communications). The rigidity of $`E_7^{}`$ can be proved easily by using a nice symmetric quadratic form of $`\mathrm{\Gamma }(𝒜^7)=E_7^{}`$ given in \[Ba94\] (see formula (2) there).

The rigid lattices $`A_1`$, $`E_6`$ and $`E_7`$ are first instances of strongly rigid lattices, i.e. such that amongst of their Delaunay polytopes there are extreme ones (see \[DGL92\]), i.e. such that any, distinct from a scaling, its affine transformation is not a Delaunay polytope. The $`1`$-simplex and unique Delaunay polytope of $`E_6`$ are only such polytopes of dimension at most $`6`$ (see \[DD01\]). $`A_1`$, $`D_4=D_4^{}`$, $`E_6`$ and $`E_7^{}`$ are rigid lattices, having unique type of Delaunay polytope, but only for $`A_1`$ and $`E_6`$ this polytope is extreme. In \[DGL92\] (see also Chapter 16 of \[DL97\]) were given 10 examples of extreme Delaunay polytopes: by one of dimension 1, 6, 7, 22. 23 and by 3 and 2 of dimension 15 and 16, respectively. We believe that there is an infinity of them. Moreover, in \[DGL92\] (see also Chapter 15 of \[DL97\]) was considered the general notion of rank of a Delaunay polytope, i.e. the number of degrees of freedom, that it has, when one deforms it affinely, so that the result remains an Delaunay polytope. Clearly, the maximal rank is (as well as maximal nrd) $`\left(\genfrac{}{}{0pt}{}{n+1}{2}\right)`$ and it is realized only by $`n`$-simplices.

In this note we describe explicitly the $`L`$-domain for the lattice $`D_n^{}`$ which is dual to the root lattice $`D_n`$. This $`L`$-domain is special and has dimension $`n`$ for odd $`n`$ and dimension 1 for even $`n`$. This special $`L`$-domain is a facet of the closure of several generic $`L`$-domains.

This work completes computation of non-rigidity degree of root lattices and their duals. In a sense, it is an addition to the work \[CS91\], where Delaunay and Voronoi polytopes of the root lattices and their duals are enumerated. We present the values of nrd for root lattices $`L`$ and their duals in the table below.

$$\begin{array}{cccccccccccc}& & & & & & & & & & & \\ L& A_1=A_1^{}& A_n& A_n^{}& D_n& D_{2m+1}^{}& D_{2m}^{}& E_6& E_6^{}& E_7& E_7^{}& E_8=E_8^{}\\ & & n2& n1& n4& m2& m2& & & & & \\ & & & & & & & & & & & \\ \text{nrd}L& 1& n+1& \frac{n(n+1)}{2}& 1& 2m+1& 1& 1& 1& 1& 1& 1\end{array}$$

## 2 The cone $`𝒢_n`$

Let $`\{e_i:iI_n\}`$ be a set of mutually orthogonal vectors of norms (i.e. of squared lengths) $`e_i^2=2\gamma _i`$, where $`I_n=\{1,2,\mathrm{},n\}`$. For $`SI_n`$, let $`e(S)=_{iS}e_i`$ and $`\gamma (S)=_{iS}\gamma _i`$. We introduce the vector $`b`$ of norm $`\alpha `$ as follows:

$$b=\frac{1}{2}\underset{iI_n}{}e_i=\frac{1}{2}e(I_n),\text{ where }b^2=\alpha =\frac{1}{2}\underset{iI_n}{}\gamma _i=\frac{1}{2}\gamma (I_n).$$
(1)
Let $`\overline{\gamma }`$ be the vector with the coordinates $`\{\gamma _i:1in\}`$. Consider the lattice $`L(\overline{\gamma })`$ generated by the vector $`b`$ and any $`n1`$ vectors $`e_i`$. If $`\gamma _i=1`$ for all $`i`$, and $`n4`$, then $`L(\overline{\gamma })=D_n^{}`$.

We take as a basis of $`L(\overline{\gamma })`$ the vector $`b`$ and the vectors $`e_i`$ for $`1in1`$. Then the coefficients of the quadratic form $`f_{\overline{\gamma }}`$ corresponding to this basis are as follows:

$$a_{ii}=e_i^2=2\gamma _i,1in1,a_{ij}=e_ie_j=0,1i,jn1,ij,$$
(2)
$$a_{nn}=b^2=\alpha ,a_{in}=e_ib=\gamma _i,1in1.$$
(3)
This form has the following explicit expression

$$f_{\overline{\gamma }}(x)=(x_nb+\underset{i=1}{\overset{n1}{}}x_ie_i)^2=\alpha x_n^2+2\underset{1}{\overset{n1}{}}\gamma _ix_i^2+2\underset{1}{\overset{n1}{}}\gamma _ix_ix_n.$$
(4)
In the basis $`\{e_i:iI_n\}`$, each vertex of $`L(\overline{\gamma })`$ has integer or half-integer coordinates.

Let $`n`$ be odd, say $`n=2m+1`$. Suppose that the parameters $`\gamma _i`$ satisfy the following $`\left(\genfrac{}{}{0pt}{}{n}{m}\right)`$ inequalities

$$\underset{iS}{}\gamma _i<\alpha ,\text{ }SI_n,\text{ }|S|=m.$$
(5)
Denote by $`𝒢_n`$ the $`n`$-dimensional domain determined in the space of variables $`\gamma _i`$, $`iI_n`$, by the inequalities (5). Since these inequalities are linear and homogeneous (recall that $`\alpha =\frac{1}{2}\gamma (I_n)`$), $`𝒢_n`$ is an open polyhedral cone. Since $`\gamma (S)+\gamma (I_nS)=2\alpha `$, the inequalities (5) imply the following inequalities

$$\gamma (T)>\alpha ,\text{ }TI_n,\text{ }|T|=m+1.$$
(6)
For a set $`T`$ of cardinality $`|T|=m+1`$, let $`T=S\{i\}`$, where $`|S|=m`$. Then (5) and (6) imply

$$\alpha <\gamma (T)=\gamma (S)+\gamma _i,\text{ i.e. }\gamma _i>\alpha \gamma (S)>0.$$

Hence the cone $`𝒢_n`$ lies in the positive orthant of $`𝐑^n`$.

Consider the closure cl$`𝒢_n`$ of the cone $`𝒢_n`$. Obviously, cl$`𝒢_n`$ is defined by the non-strict version of inequalities (5). So, using that $`2\alpha =\gamma (I_n)`$, we have

$$\mathrm{cl}𝒢_n=\{\overline{\gamma }:\gamma (S)\gamma (I_nS)0,\text{ }SI_n,|S|=m\}.$$
(7)
Note that the zero vector belongs to $`\mathrm{cl}𝒢_n`$. The automorphism group of cl$`𝒢_n`$ is isomorphic to the group of all permutations of the set $`I_n`$.

Obviously, the hyperplanes supporting facets of cl$`𝒢_n`$ are contained among the hyperplanes defined by the equalities

$$\gamma (S)=\gamma (I_nS)=\alpha =\frac{1}{2}\gamma (I_n),\text{ }SI_n,\text{ }|S|=m.$$
(8)
Note that the equality $`\gamma (S_1)=\alpha `$ can be transposed into the equality $`\gamma (S_2)=\alpha `$ by the automorphism group, for any $`S_1,S_2I_n`$ with $`|S_1|=|S_2|=m`$. Hence each of the equations of (8) determines a facet of cl$`𝒢_n`$.

###### Proposition 1

Let $`n`$ be odd and $`n=2m+15`$, i.e. $`m2`$. Then the closure of $`𝒢_n`$ has the following $`2n`$ extreme rays

$$\overline{\gamma }_q^k=\{\gamma _i=\gamma 0,iI_n\{k\},\gamma _k=2q\gamma \},\text{ }q=0,1,\text{ }kI_n.$$

Proof. Let $`\overline{\gamma }\mathrm{cl}𝒢_n`$ be fixed. Then $`\overline{\gamma }=(\gamma _1,\gamma _2,\mathrm{},\gamma _n)`$ defines a partition of the set $`I_n`$ as follows. Let the coordinates $`\gamma _i`$ take $`k`$ distinct values $`0\beta _1<\beta _2<\mathrm{}<\beta _k`$, where $`k`$ is an integer between 1 and $`n`$. For $`1jk`$, set $`S_j=\{iI_n:\gamma _i=\beta _j\}`$ and $`s_j=|S_j|`$. Then $`_{j=1}^ks_j=n=2m+1`$ and $`I_n=_{j=1}^kS_j`$ is the above mentioned partition.

Consider the values of $`\gamma (S)`$ for $`SI_n`$, $`|S|=m`$. $`\gamma (S)`$ takes a maximal value for the following sets $`S`$. Let $`j_0`$ be such that $`_{j=1}^{j_01}s_j<m+1`$, but $`_{j=1}^{j_0}s_jm+1`$. Then $`_{j=j_0+1}^ks_jm`$. Let $`S_{max}(T)=T_{j=j_0+1}^kS_j`$, where $`TS_{j_0}`$, $`|T|=t_0`$ and $`t_0:=m_{j=j_0+1}^ks_j`$. Obviously, $`\gamma (S)`$ takes the maximal value $`t_0\beta _{j_0}+_{j=j_0+1}^ks_j\beta _j`$ if $`S=S_{max}(T)`$ for any $`TS_{j_0}`$ of cardinality $`|T|=t_0`$.

For given $`\overline{\gamma }\mathrm{cl}𝒢_n`$, let $`𝒮_{max}(\overline{\gamma })`$ be the system of equations of type (8), where $`S=S_{max}(T)`$ for all $`TS_{j_0}`$ with $`|T|=t_0`$. If $`\overline{\gamma }`$ is an extreme ray of cl$`𝒢_n`$, then $`𝒮_{max}(\overline{\gamma })`$ determines uniquely up to a multiple the vector $`\overline{\gamma }`$. It is not difficult to see that $`𝒮_{max}(\overline{\gamma })`$ can uniquely determine $`\overline{\gamma }`$ only if $`k=2`$. Consider this case in detail.

If $`k=2`$, we have $`I_n=S_1S_2`$ and $`n=2m+1=s_1+s_2`$. Let $`s_2m`$, i.e. $`j_0=1`$. Then $`𝒮_{max}(\overline{\gamma })`$ consists of the following equations

$$\gamma (TS_2)=\gamma (T)+\gamma (S_2)=\alpha =\gamma (S_1T),\text{ }TS_1,\text{ }|T|=ms_2.$$

Note that $`\gamma _i`$ for $`iS_2`$ belongs to the above system only as a member of the sum $`\gamma (S_2)=_{iS_2}\gamma _i`$. Hence such a system can determine the coordinates $`\gamma _i`$, $`iS_2`$, only if $`s_2=1`$.

Now the above system implies that $`\gamma _i`$ takes the same value, say $`\gamma `$, for all $`iS_1`$. In fact, let $`i_1,i_2S_1`$, $`i_1T_1`$, $`i_1T_2`$, $`i_2T_1`$, $`i_2T_2`$, for some $`T_1,T_2S_1`$ of cardinality $`m1`$. Such $`T_1`$ and $`T_2`$ exist, since $`|T_j|=m11`$. Subtracting the equation of the above system for $`T=T_2`$ from the equation for $`T=T_1`$, we obtain the equality $`\gamma _{i_1}\gamma _{i_2}=\gamma _{i_2}\gamma _{i_1}`$, i.e. $`\gamma _{i_1}=\gamma _{i_2}`$.

In this case, the above system, where $`S_2=\{k\}`$, gives $`\gamma _k=\gamma (S_1)2\gamma (T)=s_1\gamma 2(m1)\gamma =2\gamma `$. We obtain the extreme ray $`\overline{\gamma }_1^k`$.

Now, let $`s_2>m`$, i.e. $`s_1<m+1`$ and $`j_0=2`$. A similar analysis shows that $`s_1=1`$, say $`S_1=\{k\}`$, and $`\gamma _i`$ take the same value, say $`\gamma `$, for all $`iS_2`$. This gives $`\gamma _k=0`$, and we obtain the extreme ray $`\overline{\gamma }_0^k`$. The result follows.

The facet defined by the equation $`\gamma (S)=\alpha `$, $`|S|=m`$, contains the following $`n=2m+1`$ extreme rays: $`\overline{\gamma }_0^k`$, $`kS`$, $`\overline{\gamma }_1^k`$, $`kS`$. Each facet has the following geometrical description. The $`m+1`$ rays $`\overline{\gamma }_0^k`$, $`kS`$, form an $`(m+1)`$-dimensional simplicial cone. Similarly, the $`m`$ rays $`\overline{\gamma }_1^k`$, $`kS`$, form an $`m`$-dimensional cone. Both these cones intersect by the ray $`\{\overline{\gamma }:\gamma _i=(m+1)\gamma ,iS,\gamma _i=m\gamma ,iS\}`$. Hence the cone hull of these two cones is a cone of dimension $`(m+1)+m1=2m`$. This cone is just a facet of $`𝒢_n`$ for $`n=2m+1`$.

Let $`n`$ be even, $`n=2m`$. In this case, all the inequalities (5) imply the following set of inequalities

$$\gamma (I_nS)=\gamma (T)>\alpha ,\text{ }T=I_nSI_n,|T|=m.$$

We see that this system of inequalities contradicts to the system (5). This means that the open cone $`𝒢_n`$ for even $`n=2m`$ is empty. But the solution of the set of equalities (8) is not empty. Namely, it has the solution $`\gamma _i=\gamma 0`$ for all $`iI_n`$. In other words, cl$`𝒢_n`$ is the following ray

$$\mathrm{cl}𝒢_{2m}=\{\overline{\gamma }:\gamma _i=\gamma 0,iI_{2m}\}.$$

## 3 The domain $`𝒟_n`$

Denote by $`𝒟_n`$ the domain of forms $`f_{\overline{\gamma }}`$, where $`\overline{\gamma }`$ belongs to $`𝒢_n`$.

We prove the following theorem.

###### Theorem 1

Let $`n`$ be odd, $`n=2m+1`$. The domain $`𝒟_n`$ is an $`L`$-domain. It lies in an $`n`$-dimensional space which is an intersection of $`\left(\genfrac{}{}{0pt}{}{n}{2}\right)`$ hyperplanes given by the following equalities

$$a_{ij}=0,1i<jn1,2a_{in}=a_{ii},1in1.$$
(9)
The domain $`𝒟_n`$ is cut from this space by the following inequalities

$$\underset{iS}{}a_{ii}<2a_{nn},\text{ }SI_{n1},|S|=m,$$
(10)
$$2a_{nn}<\underset{iT}{}a_{ii},\text{ }TI_{n1},|T|=m+1.$$
(11)
There is a one-to-one correspondence between $`𝒟_n`$ and the cone $`𝒢_n`$ given by the equalities (2) and (3).

In particular, the closure of $`𝒟_n`$ has $`2n`$ extreme rays $`f_0^k`$, $`f_1^k`$, $`kI_n`$, with the coefficients $`a_{ij}`$ of these forms defined as follows (where the term $`a_{kk}(f_{0,1}^k)`$ should be omitted if $`k=n`$):

$$a_{ii}(f_0^k)=2\gamma ,iI_{n1},ik,a_{kk}(f_0^k)=0,a_{nn}(f_0^k)=m\gamma ;$$

$$a_{ii}(f_1^k)=2\gamma ,iI_{n1},ik,a_{kk}(f_1^k)=4\gamma ,a_{nn}(f_1^k)=(m+1)\gamma ;$$

$`a_{ij}(f_{0,1}^k)`$ for $`ij`$ are defined by the equations (9).

The inequalities (10) and (11) define facets of the closure cl$`𝒟_n`$. All facets are domains of equivalent $`L`$-types, each having $`n`$ extreme rays $`f_1^k`$, $`kS`$, $`f_0^k`$, $`kS`$, $`SI_{n1}`$, $`|S|=m`$, or $`S=I_nT`$ and $`T`$ is as in (11).

If $`n`$ is even, $`n=2m`$, then cl$`𝒟_n`$ is one dimensional. The ray cl$`𝒟_{2m}`$ is the intersection of the $`\left(\genfrac{}{}{0pt}{}{n}{2}\right)`$ hyperplanes (9) and the $`n1`$ hyperplanes given by the following equalities

$$2a_{nn}=ma_{ii},\text{ }1in1.$$
(12)
Proof of Theorem 1 will be proceeded as follows. For a function $`f_{\overline{\gamma }}`$ given by (4), we find the Voronoi polytope. Take attention that the inequalities (10) and (11) in terms of the parameters $`\alpha `$ and $`\gamma _i`$ take the form (5) for $`nS`$ and $`nS`$, respectively. We show that the face poset of the Voronoi polytope does not change if the parameters of $`f_{\overline{\gamma }}`$ change such that they satisfy (5).

On the other hand, we show if at least one of inequalities (5) holds as equality for parameters of a function $`f_{\overline{\gamma }}`$, then the $`L`$-type of $`f_{\overline{\gamma }}`$ differs from the $`L`$-type of $`f_{\overline{\gamma }}𝒟_n`$. This will mean that $`𝒟_n`$ is an $`L`$-domain.

For to find the Voronoi polytope of $`f_{\overline{\gamma }}`$ given by (4), consider the cosets of $`2L`$ in the lattice $`L=L(\overline{\gamma })`$. Let $`v=x_nb+_{iI_{n1}}x_ie_i`$ be a vector of $`L(\overline{\gamma })`$. Then this vector belongs to the coset $`Q(S,z)`$, where $`SI_{n1}`$ is the set of indices of odd coordinates $`x_i`$ and the number $`z\{0,1\}`$ indicates the parity of the $`b`$-coordinate $`x_n`$ of the vector $`v`$. Note that the vector $`e(I_n)=2b`$ belongs to the trivial coset $`Q(\mathrm{},0)=2L`$. Hence the vectors $`e(S)`$ and $`e(I_nS)`$ belong to the same coset for any $`SI_n`$. This coset is $`Q(S,0)`$ if $`nS`$, and $`Q(I_nS,0)`$ if $`nS`$. In particular, $`e_n`$ belongs to $`Q(I_{n1},0)`$, and it is minimal in this coset. Moreover, we have $`be(S)=(be(I_nS))`$. So the $`2^n`$ vectors $`be(S)`$, $`SI_n`$, are partitioned into $`2^{n1}`$ pairs of opposite vectors.

Note that $`e^2(S)=_{iS}e_i^2=2\gamma (S)`$ and, according to (1), $`\gamma (S)+\gamma (I_nS)=2\alpha `$. Recall that $`𝒟_n`$ is the domain of $`f_{\overline{\gamma }}`$, where $`\overline{\gamma }`$ belongs to $`𝒢_n`$. Hence $`\gamma _i`$, $`iI_n`$, satisfy (5). Taking in attention (5), we see that, for $`|S|m`$, the norm of $`e(S)`$ is less than the norm of $`e(I_nS)=e(T)`$ for $`|T|m+1`$, $`S,TI_n`$.

If $`f_{\overline{\gamma }}`$ go to the boundary of $`𝒟_n`$, then the sets of minimal vectors of some cosets change. At first we describe the simple cosets which are constant on the closure of $`𝒟_n`$. The norm of minimal vectors of a coset is called also norm of the coset. These are the following cosets:

The $`n`$ cosets $`Q(\{i\},0)`$, $`iI_{n1}`$, and $`Q(I_{n1},0)`$ of norms $`2\gamma _i`$, $`iI_n`$, with minimal vectors $`e_i`$, $`iI_{n1}`$ and $`e_n=2be(I_{n1})`$, respectively.

The $`2^{n1}`$ cosets $`Q(S,1)`$ of norm $`\alpha `$ with minimal vectors $`be(S)`$, $`SI_{n1}`$.

The $`2^{n1}`$-$`n`$ non-simple cosets $`Q(S,0)`$, $`SI_n`$, $`1<|S|m`$, have norms $`\gamma (S)`$ with minimal vectors $`_{iS}\epsilon _ie_i`$, where $`\epsilon _i\{\pm 1\}`$. If $`\alpha =\gamma (S)`$, then $`|S|=m`$ and the coset $`Q(S,0)`$ contains also the vector $`_{iI_nS}\epsilon _ie_i`$.

Recall that the minimal vectors of simple cosets determine facets of the Voronoi polytope. Consider a point $`x𝐑^n`$ in the basis $`\{e_i:iI_n\}`$, $`x=_{iI_n}x_ie_i`$. Then $`x`$ belongs to the Voronoi polytope $`P`$ of $`L(\overline{\gamma })`$ if the inequalities

$$\frac{v^2}{2}xv\frac{v^2}{2}$$

hold for all minimal vectors $`v`$ of simple cosets of $`L(\overline{\gamma })`$. Using (2), (3) and the identity $`b=\frac{1}{2}_{iI_n}e_i`$, we obtain the following system of inequalities describing the Voronoi polytope of $`L(\overline{\gamma })`$:

$$\frac{1}{2}x_i\frac{1}{2},iI_n,$$
(13)
$$\frac{1}{2}\alpha \underset{iI_n}{}\gamma _i\epsilon _ix_i\frac{1}{2}\alpha ,\text{ }\epsilon _i\{\pm 1\},iI_n.$$
(14)
Here the inequality (14) is given by minimal vectors of $`Q(S,1)`$ such that $`\epsilon _i=1`$ if $`iS`$, and $`\epsilon _i=1`$ if $`iS`$.

Note that $`_{iI_n}\gamma _i\epsilon _ix_i`$ is the linear function on $`\epsilon _ix_i`$ taking maximal value if $`\epsilon _ix_i0`$ for all $`iI_n`$. Hence the right hand inequality in (14) holds as equality for a vertex $`x`$ only if $`\epsilon _ix_i>0`$ for $`x_i0`$.

An analysis of the system (13), (14) shows, that for each vertex $`x`$ there is the opposite vertex $`x`$, and $`x`$ has the following coordinates

$$x_i=\frac{1}{2}\epsilon _i,iSI_n,|S|=m,x_k=\frac{\epsilon _k}{2\gamma _k}(\alpha \gamma (S)),x_l=0,\text{ for }lI_n(S\{k\}).$$
(15)
There are $`\left(\genfrac{}{}{0pt}{}{n}{m}\right)`$ positive vertices of this type. Taking in attention signs, we obtain $`2^{m+1}\left(\genfrac{}{}{0pt}{}{n}{m}\right)`$ vertices of the Voronoi polytope. Denote the vertex (15) by $`x(k;S)`$.

The form of the vertex $`x(k;S)`$ shows that some vertices can be glued if and only if the equality $`\alpha =\gamma (S)`$ holds for some set $`S`$. If $`\alpha =\gamma (S)`$, then $`x_k(k;S)=0`$, and the $`m+1`$ vertices $`x(l;S)`$, $`lI_nS`$, are glued into one vertex. This means that if $`\alpha =\gamma (S)`$ for some $`S`$, then $`L`$-type of $`f_{\overline{\gamma }}`$ changes. So, we proved that the inequalities (5), i.e. the inequalities (10) and (11) hold for $`f𝒟_n`$.

Now, we show that $`𝒟_n`$ lies in the intersection of the hyperplanes (9). It is proved in \[BG01\] that the equations of the hyperplanes in the intersection of which an $`L`$-type domain lies are given by some linear forms on norms of minimal vectors of cosets of $`2L`$ in $`L`$. Some of such linear forms are obtained by equating norms of minimal vectors of a non-simple coset. There are $`L`$-type domains for which linear forms of last type are sufficient for to describe the space, where this $`L`$-type domain lies. This is so in our case.

In fact, it is sufficient to consider the non-simple cosets $`Q(S,0)`$ for $`|S|=2`$, i.e. to equate the norms of vectors $`e_i+e_j`$ and $`e_ie_j`$, $`i,jI_n`$. The equality $`(e_i+e_j)^2=(e_ie_j)^2`$ implies $`e_ie_j=0`$, i.e. $`a_{ij}=0`$, $`0i<jn1`$. We obtain the first equalities in (9). For $`j=n`$, we have $`e_ie_n=e_i(2be(I_{n1}))=0`$. Since $`e_ie_j=0`$, this equality is equivalent to $`2be_i=e_i^2`$. We obtain the second equalities in (9). If we set $`be_i=\gamma _i`$, $`b^2=\alpha `$, we obtain the original function $`f_{\overline{\gamma }}`$.

Now let $`n=2m`$ be even. In this case, for $`f\mathrm{cl}𝒟_n`$, all cosets of $`2L`$ in $`L`$ (excluding the cosets $`Q(S,0)`$ for $`|S|=m`$) are the same as in the odd case. But, for $`|S|=m`$, $`Q(S,0)`$ contains beside the vectors $`_{iS}\epsilon _ie_i`$ also the vectors $`_{iI_nS}\epsilon _ie_i`$. Since norms of these vectors are $`\gamma (S)`$ and $`\gamma (I_nS)`$, respectively, the equating of these norms gives the system (8). This system has the unique solution $`\gamma _i=\gamma `$ for all $`iI_n`$.

Hence, for $`n=2m`$, $`\alpha =\frac{1}{2}\gamma (I_n)=\frac{1}{2}(\gamma (S)+\gamma (I_nS))=\gamma (S)=m\gamma `$. Taking in attention (2) and (3), we can rewrite this equality as $`2a_{nn}=ma_{ii}`$ for any $`iI_n`$. So, we obtain (12). This means that cl$`𝒟_n`$ is a ray, which lies in the intersection of the hyperplanes given by the equations (10), (11) and (12). Any $`f\mathrm{cl}𝒟_n`$ is a rigid (i.e. edge) form.

So, Theorem 1 is proved.

Recall that $`L(\overline{\gamma })=D_n^{}`$ if $`\gamma _i=1`$ for all $`iI_n`$. Since, for this $`\overline{\gamma }`$, the parameters of $`f_\gamma `$ satisfy (5), where $`\alpha =m+\frac{1}{2}`$, this implies the following

Corollary $`𝒟_n`$ is the $`L`$-domain of $`D_n^{}`$. In particular, the lattice $`D_{2m}^{}`$, $`m2`$, is rigid.

Remark. Note that, for $`n`$ odd, $`n=2m+1`$, the extreme rays $`f_0^k`$ have rank $`n1=2m`$. These forms are forms of lattices isomorphic to $`\gamma D_{n1}^{}=\gamma D_{2m}^{}`$.

We saw in the proof of Theorem 1 that the Voronoi polytope of the lattice $`D_n^{}`$ is an $`n`$-cube whose vertices are cut by hyperplanes (14). A description of the Voronoi polytope of $`D_n^{}`$ can be found in \[CS91\].

Theorem 1 is a generalization of the result of \[EG01\], where the $`L`$-domain of the lattice $`D_5^{}`$ is described in detail.


text:
# Photospheric Abundances of the Hot Stars in NGC 1399 and Limits on the Fornax Cluster Cooling Flow

## 1. INTRODUCTION

The spectra of elliptical galaxies and spiral galaxy bulges exhibit a strong upturn shortward of 2700 Å, dubbed the “UV upturn.” The phenomenon was among the first major discoveries in UV extragalactic astronomy (Code 1969), and the implied existence of a hot stellar component in ellipticals contradicted the traditional picture of early-type galaxies as old, cool, passively-evolving populations. Early debate covered many possible candidates for the UV upturn origin (see O’Connell 1999), but today we know from UV spectroscopy (Ferguson et al. 1991; Brown et al. 1997) and imaging (Brown et al. 2000) that this UV emission comes from a minority population of extreme horizontal branch (EHB) stars and their progeny. These stars appear to be analogs of the subdwarf B (sdB) and O (sdO) stars of the local Galactic field, but the populations may differ considerably in the two cases.

The UV upturn is not completely understood, however. Long before the observational proof of an EHB origin, EHB stars were a strong theoretical candidate, and there have been two schools of thought regarding their nature. One possibility is that the EHB stars reside in the metal-poor tail of a wide metallicity distribution (e.g., Park & Lee 1997). Another possibility is that the EHB stars lie at high metallicity, reflecting the composition found in the cooler, dominant populations of elliptical galaxies (e.g., Brocato et al. 1990; Bressan, Chiosi, & Fagotto 1994; Greggio & Renzini 1990; Dorman, O’Connell, & Rood 1995; Horch, Demarque, & Pinsonneault 1992; Brown et al. 1997). The metallicity controversy is tied to the “second parameter” debate, which seeks to understand the role of parameters (besides metallicity) that govern HB morphology in globular clusters; possible candidates include age, mass loss, He abundance, rotation, and cluster dynamics (see Fusi Pecci & Bellazzini 1997).

Characterized by the $`m_{1550}V`$ color, the UV upturn shows surprisingly strong variation (ranging from 2.05–4.50 mag) in nearby quiescent early-type galaxies (Bertola, Capaccioli, & Oke 1982; Burstein et al. 1988), even though the spectra of ellipticals at longer wavelengths are qualitatively very similar. The $`m_{1550}V`$ color is positively correlated with the strength of Mg2 line absorption in the $`V`$ band (i.e., bluer at higher line strengths), opposite to the behavior of optical color indices (Burstein et al. 1988). Both sides in the EHB metallicity debate see this correlation as further proof of their own scenario. The metal-poor school reasons that galaxies with strong UV emission are more massive and formed earlier; such galaxies have higher mean metallicity, as tracked through optical indices, even though the upturn comes from stars in the metal-poor tail. In this scenario, age is the dominant “second parameter” driving HB morphology, with the metal-poor HB becoming more blue with age, thus producing more UV emission in massive old galaxies. In contrast, the metal-rich school argues that the trend for the HB distribution to become redder at increasing metallicity can be reversed at high metallicity, due to increased He abundance and possibly a higher mass-loss rate on the red giant branch (RGB), both of which lead to more metal-rich EHB stars and a subsequent increase in UV emission from the population. These hypotheses imply different ages for the stellar populations in these galaxies. Ages exceeding those of Galactic globular clusters are required under the metal-poor Park & Lee (1997) hypothesis, while ages as low as 8 Gyr are allowed in the metal-rich Bressan et al. (1994) model.

With this debate in mind, we obtained observations of NGC 1399 with the Far Ultraviolet Spectroscopic Explorer (FUSE; Moos et al. 2000). NGC 1399 is the giant elliptical galaxy at the center of the Fornax cluster. It has the strongest UV upturn of any quiescent elliptical galaxy measured to date ($`m_{1550}V=2.05`$ mag; Burstein et al. 1988). Although the UV upturn correlates well with optical metallicity indices, such indices track the composition of the cool population (RGB and main sequence stars). Our far-UV spectroscopy is intended to measure the photospheric abundances in the hot stars producing the UV upturn – measurements not possible in earlier spectroscopy of elliptical galaxies (Brown et al. 1997; Ferguson et al. 1991), due to the low signal-to-noise and resolution. In addition, X-ray observations suggest NGC 1399 is the host to a significant cooling flow (Bertin & Toniazzo 1995; Rangarajan et al. 1995), making it an interesting target for O VI emission measurements. We present here an analysis of our FUSE data for NGC 1399; future papers will present our observations of other ellipticals, which are not yet complete.

## 2. OBSERVATIONS

We observed NGC 1399 on 12–13 Dec 2000 and 3 Oct 2001 for 28 ksec, with 20 ksec during orbital night. The spectroscopic resolution is limited by the velocity dispersion of the galaxy (331 km s-1; Burstein et al. 1988), so we scheduled observations in the default $`30\mathrm{}\times 30\mathrm{}`$ low-resolution square aperture. The data were obtained in time-tag mode, recording the arrival time and detector position of each photon.

We reprocessed the data using the CALFUSE software, version 2.0.5. Scattered light and airglow emission lines both increase during orbital day, so we rejected daytime data. We also implemented pulse-height screening appropriate for very faint targets, by raising the minimum pulse height in CALFUSE from 0 to 4, thus discarding more of the detector dark counts with little loss of source counts.

Useful data were obtained in all detector segments, each with its own wavelength bins and instrumental resolution. Because the velocity dispersion limits the useful resolution, we coadded all of the data onto a linear wavelength scale (0.25 Å bins), with propagation of statistical errors. The bins in our coadded data are much larger than the nominal FUSE wavelength bins ($`0.01`$ Å), but they adequately sample the resolution elements in the spectrum ($`>1`$ Å for velocity-broadened features). We masked bad pixels when coadding the data, because CALFUSE does not account for small-scale flat-field features. Note that the wavelength solution in this version of CALFUSE is improved significantly, although these corrections are unimportant in our large wavelength bins. Furthermore, although fixed pattern noise can be an issue in high-resolution FUSE data, it is negligible for our large bins and coadded detector segments.

## 3. STELLAR MODELS

Our FUSE spectra offer the first clear detection of metallic absorption lines from the hot stars responsible for the UV upturn. Previously, Brown et al. (1997) fit composite models to the spectra of elliptical galaxies observed with the Hopkins Ultraviolet Telescope (HUT), including NGC 1399. Those spectra were obtained through large apertures, with coverage from 900–1800 Å, but the low signal-to-noise and resolution hampered detection of individual absorption lines from elements heavier than H. Brown et al. (1997) integrated synthetic spectra over EHB evolutionary tracks from Dorman, Rood, & O’Connell (1993), using synthetic spectra appropriate for the HUT resolution (Brown, Ferguson, & Davidsen 1996). The models reproduce well the entire far-UV spectral energy distribution, and show that the light from elliptical galaxies is dominated by a population of EHB stars and their post-HB progeny, the AGB-Manqu$`\stackrel{´}{\mathrm{e}}`$ stars (Brown et al. 1997).

To analyze the FUSE spectra of NGC 1399, we integrated the best-fitting EHB model from Brown et al. (1997), using new high-resolution synthetic spectra calculated at representative points along the EHB evolutionary path. The track assumes a small HB envelope mass (0.016 $`M_{}`$), high metallicity (\[Fe/H\]=0.71), and high main-sequence He abundance ($`Y=0.45`$). Note that approximately half of the far-UV flux in this track comes from the EHB itself, with the rest coming from the AGB-Manqu$`\stackrel{´}{\mathrm{e}}`$ phase. We computed solar abundance atmospheres using the ATLAS9 code (Kurucz 1993), under the assumption of local thermodynamic equilibrium (LTE). Synthetic spectra were then calculated using the SYNSPEC code (Hubeny, Lanz, & Jeffery 1994), again assuming LTE, but with varying abundances of C, N, Si, Fe, and Ni (holding other elements at solar abundance). Note that this code has been used previously to model FUSE and HUT data of sdB stars (Ohl, Chayer, & Moos 2000; Brown et al. 1996). We used the Kurucz (1993) line list in our spectral synthesis, including lines predicted from atomic physics but not yet measured in the laboratory. Weak lines that are experimentally unconfirmed are potentially prone to position and wavelength errors, but they provide the necessary pseudocontinuum.

No Galactic foreground reddening is applied to the models. Both the maps of Burstein & Heiles (1984) and Schlegel, Finkbeiner, & Davis (1998) suggest that $`E(BV)`$ is nearly 0 mag along a line of sight toward NGC 1399.

To compare the composite models to the FUSE data, we defined spectral indices of C, N, and Si (Table 1). Each index measures the ratio of the mean flux in a line region to that in continuum regions, selected to be free of strong features from other elements. The C and Si indices come from the ratio in one section of the spectrum, whereas the N index comes from the average in two sections.

Table 1: Composition Indices

|  | Cont. | Line | Dominant Transitions |
| --- | --- | --- | --- |
| Ion | (Å) | (Å) | (Å) |
| C III | 1170–1173 | 1174–1177 | 1174.93 1175.26 1175.59 |
|  | 1177–1180 |  | 1175.71 1175.99 1176.37 |
| N III | 986–988 | 989–993 | 989.80 991.51 991.58 |
|  | 994–996 |  |  |
| N II | 1074–1078 | 1083–1087 | 1083.99 1084.56 1084.58 |
|  | 1092–1096 |  | 1085.53 1085.55 1085.70 |
| Si III | 1105–1107 | 1107–1114 | 1108.36 1109.94 1109.97 |
|  | 1114–1116 |  | 1113.17 1113.20 1113.23 |

## 4. RESULTS

### 4.1. Photospheric Abundances

By varying the abundances in the synthetic spectra, we find that the composition that best reproduces the line indices measured in the data is C at $`0.02\pm 0.01`$ solar, N at $`0.45\pm 0.15`$ solar, and Si at $`0.13\pm 0.06`$ solar (Figure 1). The errors represent statistical uncertainties in each line index, transformed to abundance uncertainties. This result does not significantly depend upon our chosen evolutionary track; although the best-fitting track from Brown et al. (1997) has enhanced He and metallicity, our derived abundances remain unchanged (within 1$`\sigma `$) if we use other unenhanced tracks that reproduce the far-UV spectrum.

Fe and Ni lines blanket the entire FUSE region, but

Fig. 1– Sections of the FUSE NGC 1399 spectrum (black histogram). Photospheric (roman) and Galactic ISM (italics) features are easily distinguished by their widths. The composite model is shown with photospheric Fe and Ni at solar (thick grey curve) and 5 times solar (thin curve) abundance. Thick black lines at the top and bottom of the plot denote the regions used for spectral indices (Table 1). The model normalization changes by 10% at the break (the absolute flux calibration varies by $``$10% between detector segments).

do not significantly affect the derived abundances of the light elements. To demonstrate the systematic errors from Fe and Ni assumptions, we varied the abundances of Fe and Ni over a extremely large range: from 0.01 solar to 5.0 solar. Throughout this range, the derived abundances of C, N, and Si remain within their 1$`\sigma `$ statistical uncertainties. We show the same evolutionary model, when Fe and Ni are enhanced, in Figure 1; the increased Fe and Ni severely depressed the continuum, so the curve has been renormalized upward. Note that diffusion processes tend to deplete the photospheric abundances of the light elements in sdB stars (Lamontagne et al. 1985; Lamontagne, Wesemael, & Fontaine 1987), but Fe can be radiatively supported; e.g., Fe enhancement is thought to be the mechanism for the pulsating sdB stars (Charpinet et al. 1997), which show no photospheric Fe depletion (Heber, Reid, & Werner 2000).

The C and Si line indices rely upon lines that are not resonance transitions, and thus should not be contaminated from any interstellar medium (ISM) in NGC 1399 – i.e., the features in the FUSE spectrum should be completely photospheric. In contrast, the N index comes from series of resonance lines. Because the NGC 1399 spectrum is redshifted by 1424 km s-1, narrow absorption features from these same transitions in the Galactic ISM can be seen in the FUSE spectrum, offset by $`5`$ Å to the blue. If NGC 1399 has a significant cold ISM, the redshifted N features could have both ISM and photospheric contributions.

However, when looking at purely interstellar features in the FUSE data, we see no evidence for absorption from a cold ISM in NGC 1399. The N II and N III features in Table 1 can arise from both hot stellar photospheres and the ISM, but other resonance transitions will arise only from an ISM. These can be used to discern how much contamination of photospheric features we should expect from an ISM in NGC 1399. E.g., strong absorption from N I (1134.2 Å, 1134.4 Å, and 1135.0 Å) and Ca II (1135.5 Å and 1135.6 Å) from the Galactic ISM is present in the FUSE spectrum (narrow lines in Figure 1), but there is no corresponding strong absorption in the NGC 1399 frame (1140 Å, with a width of 1.3 Å). The same can be said of C II $`\lambda 1036`$ Å and O I $`\lambda 1039`$ Å (only Galactic features are present). Thus, the N II and N III spectral indices have little, if any, contribution from an ISM in NGC 1399, and should track photospheric abundances only.

### 4.2. Cooling Flow Emission

NGC 1399 lies at the center of the Fornax cluster. The X-ray luminosity of the galaxy implies a large cooling flow, with a mass deposition of at least 2 $`M_{}`$ yr-1 (Bertin & Toniazzo 1995; Rangarajan et al. 1995). Cooling flow models predict that this X-ray emission should be accompanied by O VI $`\lambda \lambda `$1031.9,1037.6 Å emission, radiated as the hot gas cools from $`10^7`$ K to $`10^5`$ K. Bregman, Miller, & Irwin (2001) did not find such emission in their FUSE observations of NGC 1404 (upper limit 0.3 $`M_{}`$ yr-1), an elliptical galaxy near the center of the Fornax cluster, but its X-ray luminosity suggests a smaller cooling flow deposition rate than that for NGC 1399 (Bertin & Toniazzo 1995). Thus, NGC 1399 offers a better opportunity for finding the warm gas from the Fornax cooling flow.

We see no O VI emission in our NGC 1399 spectrum (Figure 2); indeed, the observed flux at the location of the 1032 Å (rest frame) component is slightly lower than the expected flux from the stellar continuum. To set upper limits, we first fit an EHB model with the abundances determined above. The fit region was from 1035–1045 Å, and we included C II and O I absorption lines from our own Galactic ISM. Note that these strong Galactic features are not accompanied by equivalent features in the NGC 1399 frame, again implying a lack of a cold ISM. Once the best fit was achieved, these components were fixed, and then two components representing the O VI

Fig. 2– The NGC 1399 spectrum at the expected location of O VI emission (histogram with error bars). The curves show the stellar model (light grey), including Galactic ISM absorption (italic labels), and the 3$`\sigma `$ upper limit to the redshifted 0 VI emission (dark grey).

doublet were added. The wavelengths and widths were fixed at the NGC 1399 redshift and velocity dispersion, and the flux ratio was fixed at 2:1 for the 1032 and 1038 Å components. We find a 3$`\sigma `$ upper limit of f$`{}_{\lambda 1032}{}^{}=2.6\times 10^{15}`$ erg s-1 cm-2 (and thus f$`{}_{\lambda \lambda 1032,1038}{}^{}=3.9\times 10^{15}`$ erg s-1 cm-2 for the total O VI emission). As in Bregman et al. (2001), we assume that $`L_{\lambda 1032}=0.9\times 10^{39}\dot{M}`$ erg s-1 (a compromise between the relations for isochoric and isobaric cooling). Thus, the 3$`\sigma `$ upper limit to the mass deposition rate implied from the O VI emission is 0.14 $`M_{}`$ yr-1, for the $`3\times 3`$ kpc region sampled by the FUSE aperture (assuming $`mM=31.52`$ mag; Ferrarese et al. 2000). The limit becomes twice as large if the fit is restricted to the cleaner region around 1038 Å. Models of the X-ray luminosity imply a larger rate of $`0.4M_{}`$ yr-1 within the FUSE aperture (Rangarajan et al. 1985).

## 5. DISCUSSION

We have obtained the first clear detection of metallic absorption features in the UV upturn. The photospheric composition of the hot population in NGC 1399 appears to reflect that seen in the field sdB stars of the Milky Way. In the nearby Galactic field population, sdB and sdO stars often show moderate to severe depletions of the light elements, attributed to diffusion. Typically, N appears near solar abundance, while C and Si are depleted (Lamontagne et al. 1985; Lamontagne et al. 1987).

Several other processes could skew the photospheric abundances in EHB stars. Greggio & Renzini (1990) noted that hot HB stars of high metallicity would have very thin H envelopes; these could be wind-ejected at mass-loss rates of $`10^{10}M_{}`$ yr-1, resulting in surface compositions rich in He and N and poor in C and O. CN processing on the RGB is another possible effect (Sweigart & Mengel 1979). Assuming initial abundances at the solar C:N ratio and C+N nuclei conservation, the observed C and N abundances would imply primordial abundances of 0.12 solar, in the absence of diffusion.

Galactic field sdB stars tend to show depletion of the light elements from their original near-solar abundances, but little is known about the effects of diffusion on metal-poor EHB stars, because EHB stars in globular clusters are too distant to obtain spectroscopy easily. Diffusion is the largest uncertainty when associating the photospheric abundances of EHB stars with their primordial abundances. The photospheric abundance of a given element can be affected by both gravitational settling and radiative levitation, and thus the photospheric abundance can be either reduced or enhanced relative to the primordial abundances that drive the evolution of a star. Some calculations have shown that abundance enhancements via radiative levitation can be especially strong in the low-metallicity regime, because the absorption lines are not saturated (Michaud, Vauclair, & Vauclair 1983).

Taken at face value, our photospheric abundances do not resolve the metallicity debate for the origin of the UV upturn; they are neither extremely metal-poor nor extremely metal-rich. Tying these photospheric abundances to formation abundances will require further theoretical and observational work. In particular, UV observations of sdB stars in globular clusters are needed to demonstrate how the light elements are affected by diffusion processes in the metal-poor regime, showing whether gravitational settling or radiative levitation dominates. However, given the tendency for light element abundances to be diminished by diffusion, it seems unlikely that the UV emission in NGC 1399 comes from EHB stars with extremely low formation abundances (e.g., at $`2.4<[\mathrm{Fe}/\mathrm{H}]<1.0`$, as suggested by Park & Lee 1997); it seems more plausible that these stars started with high abundances that were then depleted, but at this point such conclusions are very speculative.

NGC 1399 is thought to be the center of a cooling flow in Fornax. Our upper limit to the O VI emission in the FUSE aperture is lower than expected from simple cooling flow models of the X-ray emission. Bregman et al. (2001) suggested that the lack of O VI emission in NGC 1404 could be due to an additional source of heating, but noted that NGC 1404 has no detected central radio source. NGC 1399 hosts a two-lobed radio source (Sadler, Jenkins, & Kotanyi 1989), which may help explain the lack of O VI emission.

Support for this work was provided by NASA through the FUSE Guest Investigator program. The authors gratefully acknowledge support from NASA grant NAS 5-9696 to the Catholic University of America. The authors wish to thank A. Sweigart and W. Landsman for useful insight and discussions.


text:
# Secure exchange of information by synchronization of neural networks

## Abstract

A connection between the theory of neural networks and cryptography is presented. A new phenomenon, namely synchronization of neural networks is leading to a new method of exchange of secret messages. Numerical simulations show that two artificial networks being trained by Hebbian learning rule on their mutual outputs develop an antiparallel state of their synaptic weights. The synchronized weights are used to construct an ephemeral key exchange protocol for a secure transmission of secret data. It is shown that an opponent who knows the protocol and all details of any transmission of the data has no chance to decrypt the secret message, since tracking the weights is a hard problem compared to synchronization. The complexity of the generation of the secure channel is linear with the size of the network.

(1) Department of Physics, Bar Ilan University, 52900 Ramat Gan, Israel

(2) Institut für Theoretische Physik, Universität Würzburg, Am Hubland, D-97074 Würzburg, Germany

3. 9. 2001

PACS numbers: 87.18.Sn,89.70.+c

The ability to build a secure channel is one of the most challenging fields of research in modern communication. Since the secure channel has many applications, in particular for mobile phone, satellite and internet-based communications, there is a need for fast, effective and secure transmission protocols . Here we present a novel principle of a cryptosystem based on a new phenomenon which we observe for artificial neural networks.

The goal of cryptography is to enable two partners to communicate over an insecure channel in such a way that an opponent cannot understand and decrypt the transmitted message. In a general scenario, the message is encrypted by the sender through a key $`E_k`$ and the result, the ciphertext, is sent over the channel. A third party, eavesdropping on the channel, should not determine what the message was. However, the recipient who knows the encryption key can decrypt the ciphertext using his private key $`D_k`$.

In a private–key system the recipient has to agree with the sender on a secret key $`E_k`$, which requires a hidden communication prior to the transmission of any message. In a public–key system, on the other side, the key $`E_k`$ is published and a hidden communication is not necessary. Nevertheless, an opponent cannot decrypt the transmitted message since it is computationally infeasible to invert the encryption function without knowing the key $`D_k`$. In a key-exchange protocol, both partner start with private keys and transmit – using a public protocol – their encrypted private keys which, after some transformations, leads to a common secret key. In most applications a public-key system is used which is based on number theory where the keys are long integers .

In this report we suggest a novel cryptosystem. It is a key-exchange protocol which does neither use number theory nor a public key, but it is based on a learning process of neural networks: The two participants start from a secret set of vectors $`E_k(0)`$ and $`D_k(0)`$ without knowing the key of their partner. By exchanging public information the two keys develop to a common time dependent key $`E_k(t)=D_k(t)`$, which is used to encrypt and decrypt a given message. An opponent who knows the algorithm and observes any exchange of information is not able to find the keys $`E_k(t)`$ and $`D_k(t)`$. Our method is based on a new phenomenon presented here: Synchronization of neural networks by mutual learning .

Simple models of neural networks describe a wide variety of phenomena in neurobiology and information theory . Artificial neural networks are systems of elements interacting by adaptive couplings which are trained from a set of examples. After training they function as content addressable associative memory, as classifiers or as prediction algorithms.

In this report we present a new phenomenon: Two feedforward networks can synchronize their synaptic weights by exchanging and learning their mutual outputs for given common inputs. Surprisingly, synchronization is fast; the number of bits required to achieve perfect alignment of the weights is lower than the number of components of the weights. After synchronization, the synaptic weights define the common time dependent private key $`E_k(t)=D_k(t)`$. With respect to possible applications we find that first, tracking the weights of one of the networks by the opponent is a hard problem. Although we were not able to find a mathematical proof, our simulations, in addition to arguments based on analytic results on neural networks, give clear evidence that our key exchange protocol is secure . Second, the complexity of our cryptosystems scales linearly with the size of the network (=number of bits of the keys). In summary, from this new biological mechanism one can construct efficient encryption systems using keys which change permanently.

This phenomenon, as well as the corresponding applications in cryptography, can be extended to a system of several partners communicating with each other, as well as to other tasks relying on a secure channel . Since synchronization is a subject of recent research in neuroscience too , we believe that our bridge between the theory of neural networks and cryptography may help to understand communication between parts of biological neuronal or genetic networks.

In the following we introduce and investigate a simple model which shows the properties sketched above. The architecture used by the recipient and the sender is a two-layered perceptron, exemplified here by a parity machine (PM) with $`K`$ hidden units. More precisely, the size of the input is $`KN`$ and its components are denoted by $`x_{kj},k=1,2,\mathrm{},K`$ and $`j=1,\mathrm{},N`$. For simplicity, each input unit takes binary values, $`x_{kj}=\pm 1`$. The $`K`$ binary hidden units are denoted by $`y_1,y_2,\mathrm{},y_K`$. Our architecture is characterized by non-overlapping receptive fields (a tree), where the weight from the j$`th`$ input unit to the k$`th`$ hidden unit is denoted by $`w_{kj}`$, and the output bit $`O`$ is the product of the state of the hidden units (see Fig. 1). For simplicity we discuss PMs with three hidden units $`K=3`$. We use integer weights bounded by $`L`$, i.e. $`w_{kj}`$ can take the values $`L,L+1,\mathrm{},L`$.

The secrete information of each of the two partners are the initial values for the weights, $`w_{kj}^S`$ and $`w_{kj}^R`$, for the sender and the recipient, respectively. It consists of $`6N`$ integer numbers, $`3N`$ of the recipient and $`3N`$ of the sender. Sender and recipient do not know the initial numbers of their partners, which are used to construct the common secret key.

Now each network is trained with the output of its partner. At each training step, for the synchronization as well as for the encryption/decryption step, a new common public input vector $`(x_{kj})`$ is needed for both the sender and the recipient. For a given input, the output is calculated in the following two steps. In the first one, the state of the three hidden units, $`y_k^{S/R},k=1,2,3`$, of the sender and the recipient are determined from the corresponding fields

$$y_k^{S/R}=\text{sign}[\underset{j=1}{\overset{N}{}}w_{kj}^{S/R}x_{kj}]$$
(1)
In the case of zero field, $`w_{kj}^{S/R}x_{kj}=0`$, the sender/recipient sets the hidden unit to $`1/1`$. In the next step the output $`O^{S/R}`$ is determined by the product of the hidden units, $`O^{S/R}=y_1^{S/R}y_2^{S/R}y_3^{S/R}`$.

The sender is sending its output (one bit) to the recipient, the recipient is sending its output to the sender and both networks are trained with the output of its partner. In case that they do not agree on the current output, $`O^SO^R<0`$, the weights of the sender/recipient are updated according to the following Hebbian learning rule .

$`\text{if}\left(O^{S/R}y_k^{S/R}>0\right)`$ then $`w_{kj}^{S/R}=w_{kj}^{S/R}O^{S/R}x_{kj}`$
$`\text{if}\left(|w_{kj}^{S/R}|>L\right)`$ then $`w_{kj}^{S/R}=\text{sign}(w_{kj}^{S/R})L`$ (2)
Only weights belonging to the one (or three) hidden units which are in the same state as that of their output unit are updated, in each one of the two networks. Note that by using this dynamical rule, the sender is trying to imitate the response of the recipient and the recipient is trying to imitate the one of the sender.

There are three main ingredients in our model which are essential for a secure key exchange protocol: First, from the knowledge of the output, the internal representation of the hidden units is not uniquely determined because there is a four fold degeneracy (for the output $`+1`$ there are four internal representations for the three hidden units $`(1,1,1),(1,1,1),(1,1,1),(1,1,1)`$). As a consequence, an observer cannot know which of the weight vectors is updated according to equation (2). Second, we have chosen the parity machine since in the case of static weights, it is known that an opponent cannot obtain any knowledge about the rule if he is trained with less than $`\alpha (L)N`$ random examples (where $`\alpha (L)=2.63`$ for $`L=3`$ and for more details see ). This analytic result favours the PM over other multilayer networks. Third, since each component is bounded by $`L`$, an observer cannot invert the sum of equation (2); the network forgets . As a consequence of these three ingredients, the initial weight vectors cannot be recovered from the knowledge of the time dependent synchronized keys. All three of these mechanisms – hidden units, PM as well as bounded weights – make the problem hard for any observer.

We find that the two PMs learning from each other are able to synchronize, at least for some parameters $`K,L`$ and $`N`$ . Our simulations show that after an initial relatively short transient time the two partners align themselves into antiparallel states. It is easy to verify from our learning rule that as soon as the two networks are synchronized they stay so forever. The number of time steps to reach this state depends on the initial weight vectors and on the sequence of random inputs, hence it is distributed. Fig. 2 shows the distribution of synchronization time obtained from at least $`1000`$ samples. It is evident that two communicating networks synchronize in a rather short time. The average synchronization time $`t_{av}`$ decreases with increasing size $`N`$ of the system, see Fig. 3; it seems to converge to $`t_{av}410`$ for infinitely large networks. Surprisingly, in the limit of large $`N`$ one needs to exchange only about $`400`$ bits to obtain agreement between $`3N`$ components. However, one should keep in mind that the two partners do not learn the initial weights of each other, they just are attracted to a dynamical state with opposite weight vectors.

As soon as the weights of the sender and the recipient are antiparallel the public initialization of our private-key cryptosystem is terminated successfully and the encryption of the message starts. Now there are two possibilities to choose an algorithm: First, use a conventional encryption algorithm, for example a stream cipher like the well-known Blum-Blum-Shub bit generator . In this case the seed for this pseudo-random number generator is constructed from our weight vector after synchronization. Second, use the PM itself for a stream cipher by multiplying its output bit with the corresponding bit of the secret data.

In the case of the PM, the complexity of the encryption/decryption processes scales linearly with the size of the transmitted message, whereas the complexity of the synchronization process does not scale with the size of the network. Hence our construction is a linear cryptosystem .

Now we examine a possible attack on our cryptosystem. The opponent eavesdropping on the line knows the algorithm as well as the actual mutual outputs, hence he knows in which time steps the weights are changed. In addition, the opponent knows the input $`x_{kj}`$ as well. However, the opponent does not know the initial conditions of the weights of the sender and the recipient. As a consequence, even for the synchronized state, the internal representations of the hidden units of the sender and the recipient are hidden from the opponent and he does not know which are the weights participating in the learning step. For random inputs all four internal representations appear with equal probability in any stage of the dynamical process, hence for $`t`$ training steps there are $`4^t`$ possibilities to select internal representations.

Therefore, on the time scale of synchronization the observing network has no chance to obtain complete knowledge about the other two networks. We have simulated attacks of an observer, assuming that the most effective algorithm is a network which has identical architecture to the recipient: A PM with the same learning rule and parameters as described above. The observing network is trained with the input vector and output bit of the sender, and the training step (2) is performed only if sender and recipient disagree with each other. Note that the sender does not react to the output of the opponent, which results in a large noise to signal ratio compared with the recipient .

The learning rule (2) may be considered for each component of the weight vectors as a kind of biased random walk with reflecting boundaries. Therefore, for very long times, the observer may take the weight vector of the other network by chance. The distribution of the ratio between the time $`t_{sync}`$ the sender and recipient need to synchronize and the learning time $`t_{learn}`$ the opponent needs for complete overlap is shown in Fig. 4. For $`N=101`$ the average learning time is a factor of about $`125`$ larger than the corresponding synchronization time. In addition, with increasing system size the tail of the distribution for larger ratios is reduced.

Hence the time to synchronize by chance is very long and in the example discussed here it is of order $`O(10^5)`$ . The heart of our cryptosystem is that synchronization is a much simpler task than tracking by an observer.

This principle is also supported by the following observation. Assume that the observer has already some knowledge about one of the networks, he knows $`M`$ out of $`N`$ components for each of the three hidden units. We have measured whether the observer succeeds to synchronize within 2000 time steps – a time much longer than the average synchronization time for the two communicating networks. For $`M=N`$ the observer is parallel to one of the networks and remains so forever. But already for $`M=N1`$ we find a high probability that the weight vector of the observer separates completely from its almost parallel alignment. For smaller values of $`M`$ this probability decreases fast to zero. Surprisingly, even in case the observer has almost complete knowledge about the two partners, he does not succeed to achieve complete information from learning examples. This fact reduces the probability of the opponent to imitate one of the communicating networks using an ensemble of PMs.

Our key exchange protocol can be generalized to include Bit-Packages as is briefly described below. An important issue for the implementation of our cryptosystem is to accelerate the synchronization process from hundreds of time steps to a few dozens while keeping the security of our channel. Surprisingly, both of these two goals can be achieved simultaneously sending bit packages (BP). In this scenario the process contains the following steps: (a) The sender and the recipient generate $`B>1`$ common inputs. (b) The sender and the recipient calculate the output of their PM for the set of $`B`$ inputs and store the $`B`$ sets of the corresponding values $`y_{ki}(i=1,\mathrm{},B)`$ of the hidden units (the internal representations) (c) The transmission of mutual information; the sender/recipient sends a package consisting of $`B`$ bits ($`b_i^{S/R}`$) to the recipient/sender. (d) The sender and the recipient are updating their weights using the same learning rule as for $`B=1`$: In case that bit $`b_i^Sb_i^R`$ the learning process is taking place as before using the corresponding internal representations. The synchronization time is dramatically reduced, as is shown in Fig. 5. For instance, for $`N=21`$, $`K=3`$, $`L=3`$, synchronization is achieved after 12 bit packages if the size of the package is larger than $`B32`$.

Other extensions of our method as well as the analytical calculation of the distribution of $`t_{sync}`$, $`t_{learn}`$ for various $`K,L`$ and continuous weights, and version space of the PMs which are consistent with the training set will be dicussed in .

Finally, we want to remark that synchronization is a subject of recent research in neuroscience, where for instance, in experiments on cats and monkeys one has found that the spike activity of neurons in the visual cortex has correlations which depend on the kind of optical stimulus shown to the animal . The phenomenon described here suggests that synchronization can be used by biological neuronal networks or by networks of the immune system to exchange secure information between different parts of an organism.