text:
# 1 Introduction

## 1 Introduction

In a recent paper , Verlinde argued that the Cardy formula , describing the entropy of a certain conformal field theory (CFT) in $`1+1`$ dimensions, can be generalized to any dimension, leading to the so-called Cardy-Verlinde formula. Consider a certain CFT residing in an $`(n+1)`$-dimensional spacetime with the metric

$$ds^2=dt^2+R^2d\mathrm{\Omega }_n^2,$$
(1.1)
where $`d\mathrm{\Omega }_n^2`$ denotes the line element of a unit $`n`$-dimensional sphere and $`R`$ is the radius of the sphere. It was proposed that the entropy of the CFT in the spacetime (1.1) can be related to its total energy $`E`$ and Casimir energy $`E_c`$ as

$$S=\frac{2\pi R}{\sqrt{ab}}\sqrt{E_c(2EE_c)}.$$
(1.2)
Here $`a`$ and $`b`$ are two positive parameters. For strongly coupled CFTs with the AdS duals, which implies that the CFTs are in the regime of supergravity duals, $`ab`$ is fixed to be $`n^2`$ exactly. Thus one obtains the Cardy-Verlinde formula

$$S=\frac{2\pi R}{n}\sqrt{E_c(2EE_c)}.$$
(1.3)
Indeed, this formula holds for various kinds of asymptotically AdS spacetimes whose boundary metric is of the form (1.1): AdS Schwarzschild black holes ; charged AdS black holes  and Taub-Bolt AdS spacetimes , whose thermodynamics corresponds to that of different CFTs . The Cardy-Verlinde formula (1.3) holds even for the AdS Kerr black holes , whose boundary is a rotating Einstein universe.111For more references discussing the Cardy-Verlinde formula see -.

More recently, it has been proposed that defined in a manner analogous to the AdS/CFT correspondence, quantum gravity in a de Sitter (dS) space is dual to a certain Euclidean CFT living on a spacelike boundary of the dS space  (see also earlier works -). Following the proposal, some investigations on the dS space have been carried out recently -.

It is well known that in de Sitter space, there is a cosmological horizon, which has the similar thermodynamic properties like the black hole horizon . According to the dS/CFT correspondence, it might be expected that as the case of AdS black holes , the thermodynamics of cosmological horizon in asymptotically dS spaces can be identified with that of a certain Euclidean CFT residing on a spacelike boundary of the asymptotically dS spaces. Thus it is of great interest to see whether the entropy of the cosmological horizon can be described by the Cardy-Verlinde formula (1.3). This is the purpose of the present paper.222see  for another attempt in this direction. The dynamics of a brance in a Schwarzschild-dS spacetime has been discussed in .

To this end, we have to first face a difficulty to calculate some conserved charges including the mass (total energy) of gravitational field of the asymptotically dS spaces. In the spirit of the dS/CFT correspondence, these conserved charges of gravitational field can be identified with those of corresponding Euclidean CFT. The difficulty arises due to the absence of the spatial infinity and the globally timelike Killing vector in an asymptotically dS space. In this paper we will follow the prescription recently proposed in  to calculate the mass of gravitational field of asymptotically dS spaces (and then the energy of corresponding CFTs). It is found that indeed the entropy of cosmological horizon in some asymptotically dS spaces can be described in terms of the Cardy-Verlinde formula.

The organization of the paper is as follows. In the next section we first discuss the case of Schwarzschild-dS solutions. In Sec. 3 we will consider the case of topological dS solutions presented in . We conclude with a discussion in Sec. 4.

## 2 Schwarzschild-dS Solutions

Consider an $`n+2`$-dimensional Schwarzschild-dS spacetime333In three dimensions, it has been already shown that the entropy of Schwarzschild (Kerr)-de Sitter solution can be expressed by Cardy formula, for example see  and references therein. So in this paper we consider the case of spacetime dimension $`n+24`$.

$$ds^2=f(r)dt^2+f(r)^1dr^2+r^2d\mathrm{\Omega }_n^2,$$
(2.1)
where

$$f(r)=1\frac{\omega _nm}{r^{n1}}\frac{r^2}{l^2},\omega _n=\frac{16\pi G}{nVol(S^n)}.$$
(2.2)
Here $`m`$ is an integration constant, $`Vol(S^n)`$ denotes the volume of a unit $`n`$-sphere $`d\mathrm{\Omega }_n^2`$, and $`G`$ is the gravitational constant. The Schwarzschild-dS spacetime (2.1) is a solution of Einstein equations with a cosmological constant $`\mathrm{\Lambda }=n(n+1)/2l^2`$ in $`n+2`$ dimensions. When the Schwarzschild parameter $`m`$ vanishes, the solution (2.1) reduces to the dS solution, which has a cosmological horizon at $`r=l`$. The horizon has the associated Hawking temperature $`T_{\mathrm{HK}}=1/2\pi l`$ and entropy $`S=l^nVol(S^n)/4G`$ . When $`m`$ increases with $`m>0`$, a black hole horizon occurs and increases in size with $`m`$, while the cosmological horizon shrinks. Finally the black hole horizon and cosmological horizon coincide with each other when

$$m_N=\frac{2}{\omega _n(n+1)}\left(\frac{n1}{n+1}l^2\right)^{(n1)/2}.$$
(2.3)
This is the Nariai black hole, the maximal black hole in dS space.

When $`m0`$, the cosmological horizon $`r_c`$ of the solution (2.1) is determined by the maximal root of the equation $`f(r)=0`$. The associated Hawking temperature and entropy to the cosmological horizon are

$`T_{\mathrm{HK}}={\displaystyle \frac{1}{4\pi r_c}}\left((n+1){\displaystyle \frac{r_c^2}{l^2}}(n1)\right),`$
$`S={\displaystyle \frac{r_c^nVol(S^n)}{4G}}.`$ (2.4)
In Ref. , using surface counterterm method the authors calculated the mass of gravitational field of the spacetime (2.1) in four and five dimensions 444For the case in higher dimensions, more counterterms than those given in  are needed.. It was found that

$$M_4=m,M_5=\frac{3\pi l^2}{32G}m.$$
(2.5)
Thus a pure dS4 has a vanishing mass, while the Nariai black hole has the mass $`M_4=l/3\sqrt{3}G`$ in four dimensions. The pure dS5 has the mass $`3\pi l^2/32G`$ and the Nariai black hole in five dimensions has a vanishing mass. This is a surprising result: The mass of black hole solution is always less than that of pure dS space in corresponding dimensions. However, it was argued that this result is consistent with the Bousso’s observation on the asymptotically dS space and the dS/CFT correspondence . The Bousso’s observation is that the entropy of dS space is an upper bound for the entropy of any asymptotically dS spaces . Furthermore, if the dS/CFT correspondence exists, the mass of asymptotically dS spaces should be translated into the energy of a dual Euclidean CFT. Generically it is expected that such a field theory has entropy increasing with energy. The result (2.5) has precisely this property .

The nonvanishing mass of pure dS5 space is reminiscent of the nonvanishing mass of pure AdS5 in the global coordinates . The latter is shown to be just the Casimir energy of $`𝒩=4`$ SYM theory on a $`3`$-sphere. Such a Casimir energy is expected to exist for pure AdS space in odd dimensions . In the closest analogy of the case in AdS spaces, the nonvanishing mass of pure dS spaces is also expected to exist in odd dimensions in some coordinates (cf. ). These nonvanishing masses can be viewed as the Casimir energies of CFTs dual to the pure dS spaces.

Note that in checking the Cardy-Verlinde formula for the asymptotically AdS spaces, the Casimir energy of pure AdS spaces did not be included 555If one wants to include the Casimir energy of the pure AdS space, this part must be subtracted from the total energy in the Cardy-Verlinde formula, in a way as we did in the case of charged AdS black holes .. In a manner analogous to this, we will not also consider the Casimir energy corresponding to the pure dS spaces in the present context. As a result, we obtain the gravitational mass of the Schwarzschild-dS solution (2.1) from (2.5666For other higher dimensions than $`4`$ and $`5`$, we expect that the following result holds also.

$$E=m=\frac{r_c^{n1}}{\omega _n}\left(\frac{r_c^2}{l^2}1\right),$$
(2.6)
where $`r_c`$ denotes the cosmological horizon. According to the dS/CFT correspondence, the gravitational mass (2.6) is just the energy of the dual CFT. Thus, in this case the energy of CFT is negative, a unpleasant result777We have more discussions on this below.. Despite of the non-positive definiteness of the energy, we first follow naively the steps as the case of AdS black holes to check whether the entropy (2) of the cosmological horizon can be expressed or not by the Cardy-Verlinde formula (1.3).

Following Ref. , the Casimir energy $`E_c`$ can be calculated using the formula

$$E_c=n(E+pVT_{\mathrm{HK}}S),$$
(2.7)
where $`p`$ is the pressure and $`V`$ is the volume. Since we are considering a CFT, so we have $`p=E/nV`$. Corresponding to the Schwarzschild-dS solution (2.1), the dual Euclidean CFT resides on the space

$$ds^2=dt^2+l^2d\mathrm{\Omega }_n^2.$$
(2.8)
Note that here $`t`$ denotes a spacelike coordinate, rather than a timelike coordinate. Substituting (2) and (2.6) into (2.7), we obtain

$$E_c=\frac{2nr_c^{n1}Vol(S^n)}{16\pi G}.$$
(2.9)
This Casimir energy is negative, which is reminiscent of the case of the hyperbolic AdS black holes, there the Casimir energy is also negative . The negative Casimir energy indicates that the dual CFT is not unitary. This agrees with Ref. , in which it was argued that the dual Euclidean CFT to the dS space is not unitary.

Given the Casimir energy (2.9), it is easy to check that the entropy $`S`$ (2) can be expressed as follows

$$S=\frac{2\pi l}{n}\sqrt{|E_c|(2EE_c)}.$$
(2.10)
This form is completely the same as the Cardy-Verlinde formula in the case of hyperbolic AdS black holes . Note that this formula (2.10) also holds for pure dS spaces. In that case, $`E=0`$ and

$$E_c=\frac{2nl^{n1}Vol(S^n)}{16\pi G},$$
(2.11)
the expression (2.10) then gives the entropy of pure dS spaces. The result (2.10) looks fine. But there are two points to be understood. The first is that in this formula the energy of corresponding CFT is negative888If replacing $`E=m`$ by $`E=m`$ in (2.6), one then cannot arrive at the Cardy-Verlinde formula (2.10).. The other is that the formula (2.10) gives us the entropy of cosmological horizon only, and does not apply to the black hole horizon. As we stated above, when $`mm_N`$, except for the cosmological horizon, a black hole horizon $`r_+`$ occurs, which has the Hawking temperature $`\stackrel{~}{T}_{\mathrm{HK}}`$ and entropy $`\stackrel{~}{S}`$

$`\stackrel{~}{T}_{\mathrm{HK}}={\displaystyle \frac{1}{4\pi r_+}}\left((n1)(n+1){\displaystyle \frac{r_+^2}{l^2}}\right),`$
$`\stackrel{~}{S}={\displaystyle \frac{r_+^nVol(S^n)}{4G}}.`$ (2.12)
For the Schwarzschild-dS solution, the total entropy should be the sum of black hole horizon entropy $`\stackrel{~}{S}`$ and cosmological horizon entropy $`S`$. But we cannot obtain a similar formula like (2.10) for the entropy $`\stackrel{~}{S}`$ of black hole horizon. Therefore further investigation is needed for the Schwarzschild-dS spacetime. In order to avoid the two points, in the next section we consider the topological dS solutions presented in .

## 3 Topological dS Solutions

In Ref. , an interesting conjecture was put forward, which states that any asymptotically de Sitter space whose mass exceeds that of de Sitter contains a cosmological singularity. To check this conjecture, in Ref.  Myung, Zhang and the present author present a solution named topological de Sitter solution. The solution can be described by the metric

$$ds^2=f(r)dt^2+f(r)^1dr^2+r^2\gamma _{ij}dx^idx^j,$$
(3.1)
where

$$f(r)=k+\frac{\omega _nm}{r^{n1}}\frac{r^2}{l^2},\omega _n=\frac{16\pi G}{nVol(\sigma )},$$
(3.2)
$`m`$ is an integration constant which is supposed to be positive, and $`\gamma _{ij}dx^idx^j`$ denotes the line element of an $`n`$-dimensional hypersurface $`\sigma `$ with constant curvature $`n(n1)k`$ and volume $`Vol(\sigma )`$. Without loss of generality, the constant $`k`$ can be set to $`1`$, $`0`$ or $`1`$.

The solution (3.1) is asymptotically de Sitter. But because of $`m0`$, black hole horizon is not present in this solution, instead a naked singularity at $`r=0`$ occurs when $`m0`$. Although so, the solution has a cosmological horizon $`r_c`$, which is the root of the equation, $`f(r)=0`$. And the associated Hawking temperature $`T_{\mathrm{HK}}`$ and entropy $`S`$ to the cosmological horizon are

$`T_{\mathrm{HK}}={\displaystyle \frac{1}{4\pi r_c}}\left((n+1){\displaystyle \frac{r_c^2}{l^2}}(n1)k\right),`$
$`S={\displaystyle \frac{r_c^nVol(\sigma )}{4G}}.`$ (3.3)
When $`k=1`$, the solution (3.1) has the same form as that of the Schwarzschild-dS solution if $`m`$ is replaced by $`m`$ in (2.1). Due to this, following the prescription in Ref. , it is easy to show that the gravitational masses of the solutions (3.1) in four and five dimensions are

$$M_4=m,M_5=\frac{3\pi l^2}{32G}+m,$$
(3.4)
respectively. When $`k=0`$, the gravitational mass was found to be

$$M=m,$$
(3.5)
while when $`k=1`$, the gravitational masses in four and five dimensions are

$$M_4=m,M_5=\frac{3l^2Vol(\sigma )}{64\pi G}+m.$$
(3.6)
For these cases, the gravitational masses are always larger than those of de Sitter spaces ($`m=0`$) in corresponding dimensions. As a result, we verify the conjecture in Ref.  within these examples. For more details see .

Corresponding to the solution (3.1), the dual Euclidean CFT resides on the space

$$ds^2=dt^2+l^2\gamma _{ij}dx^idx^j.$$
(3.7)
Once again, here the coordinate $`t`$ is a spacelike one. As the case of asymptotically AdS spaces, we neglect the nonvanishing mass of pure dS spaces, as those in (3.4) and (3.6), or subtract the nonvanishing energy from the total energy of corresponding CFTs, the energy of the Euclidean CFTs dual to the solution (3.1) is then

$$E=m=\frac{r_c^{n1}}{\omega _n}\left(\frac{r_c^2}{l^2}k\right).$$
(3.8)
Note that here the energy of CFTs is always positive, different from the case (2.6) of the Schwarzschild-dS solution. Substituting (3) and (3.8) into (2.7), we get the Casimir energy $`E_c`$

$$E_c=\frac{2nkr_c^{n1}Vol(\sigma )}{16\pi G}.$$
(3.9)
When $`k=0`$, the Casimir energy vanishes, as the case of asymptotically AdS spaces. This is expected since the thermodynamic quantities of CFTs in a Ricci flat space are conformal invariant, there is no finite volume effect. When $`k=\pm 1`$, we see from (3.9) that the sign of the energy is just contrast to the case of asymptotically AdS space , there it was found that for a hyperbolic space with $`k=1`$, the Casimir energy is negative , while it is positive for the $`k=1`$ case .

With the Casimir energy (3.9), one can easily see that the entropy (3) of the cosmological horizon can be expressed in a form of the Cardy-Verlinde formula as

$$S=\frac{2\pi l}{n}\sqrt{|E_c|(2EE_c)},$$
(3.10)
when $`k=\pm 1`$. To accommodate the case of $`k=0`$, we can rewrite (3.10) as

$$S=\frac{2\pi l}{n}\sqrt{|E_c/k|(2EE_c)}.$$
(3.11)
As a result, we show that indeed the entropy of the cosmological horizon in the topological dS solution (3.1) can be described by the Cardy-Verlinde formula (1.3). This also further provides evidence in favour of the dS/CFT correspondence.

## 4 Conclusions

In the dS/CFT correspondence, we have investigated the question of whether the entropy of cosmological horizon in asymptotically dS spaces can be described by the Cardy-Verlinde formula, which was established in the AdS/CFT correspondence . For the Schwarzschild-dS solution, although the gravitational mass, calculated in the prescription of Ref. , of the solution (and then the energy of the dual Euclidean CFT) is always negative, we have found that the entropy of the cosmological horizon (and then the entropy of the dual CFT) indeed can be expressed in terms of a form \[see (2.10)\] of the Cardy-Verlinde formula. However, we cannot find a similar formula for the entropy of black hole horizon in the Schwarzschild-de Sitter spacetime. Therefore further study is needed for the Schwarzschild-dS solution. When the Schwarzschild parameter $`m`$ vanishes, our result (2.10) precisely reproduces the entropy of pure dS spaces.

In the topological dS solutions , which have a cosmological horizon and a naked singularity at $`r=0`$. The gravitational mass is always positive. The entropy associated to the cosmological horizon was found to obey the Cardy-Verlinde formula \[see (3.11)\]. This result is in favour of the dS/CFT correspondence. One point which should be further investigated is whether the Euclidean CFT dual to the topological dS solution can describe the naked singularity in this asymptotically dS spacetime. In other words, whether is there a well-defined Euclidean CFT dual to an asymmetrically dS space which contains a naked singularity?

It would be of interest to extend the discussions made in this paper to the case of asymptotically dS spaces with rotation and/or charges.

## Acknowledgments

The author thanks Y.S. Myung for useful discussions. This work was supported in part by a grant from Chinese Academy of Sciences.


format: text

text:
# On the effective temperature scale of O stars

## 1 Introduction

As a significant fraction of the flux of O stars is emitted in the inaccessible Lyman continuum ($`\lambda <912`$ Å) reliable direct determinations of their effective temperatures are not possible. Indirect methods, primarily based on atmospheric modeling, are therefore employed (e.g. Böhm-Vitense bv81 (1981), Crowther pac98 (1997)). Given the need for a detailed treatment of non-LTE effects and the presence of stellar winds (Kudritzki & Hummer kh90 (1990)), a complete modeling of such atmospheres including also the effects of numerous metal-lines (“line blanketing”) remains a complex task (cf. Schaerer & Schmutz ss94 (1994), Hillier & Miller hm98 (1998), Pauldrach et al. p01 (2001)).

For these reasons, most published spectral analysis have so far been based on simple non-LTE models. For example, the most recent calibration of stellar parameters of O and early B type stars of Vacca et al. (vacca (1996), hereafter VGS96) is based only on results from plane parallel, pure Hydrogen and Helium (H-He) non-LTE models. Their derived temperature scale for O stars is found to be significantly hotter than most earlier calibrations (see references in VGS96). Such differences lead to non-negligible changes in the fundamental parameters of O stars — e.g. luminosities, Lyman continuum fluxes etc. — when estimated from spectral types. Accurate calibrations are crucial for various astrophysical topics, such as comparisons with stellar evolution models, determinations of the initial mass function and cluster ages, studies of H ii regions, and others.

Indications for a decrease of $`T_{\mathrm{eff}}`$ due to line blanketing effects have been found since the first non-LTE $`+`$ wind modeling attempts by Abbott & Hummer (ah85 (1985), and subsequent investigations based on the same “wind blanketed” models), the improved models of Schaerer & Schmutz (ss94 (1994)) and Schmutz (s98 (1998)), and the fully-blanketed plane parallel non-LTE models of Hubeny et al. (hub98 (1998)). Similar indications are obtained by Fullerton et al. (ful00 (2000)) from recent modeling of FUSE spectra with the code of Pauldrach et al. (p01 (2001)) and by Crowther et al. (cro01 (2001)).

The effective temperature scale of O stars is revised here based on the recent CMFGEN code of Hillier & Miller (hm98 (1998)), which treats the problem of a non-LTE line blanketed atmosphere with a stellar wind in a direct way, thereby avoiding possible shortcomings due to opacity sampling techniques employed by Schaerer & Schmutz (ss94 (1994)), Schmutz (s98 (1998)), and Pauldrach et al. (p01 (2001)). First results on the dwarf sequence are presented here. A more detailed account including all luminosity classes will be presented in a subsequent publication.

In Sect. 2 we describe our method and the calculated models. The results, their dependence on model/stellar parameters, and first comparison with observations are presented in Sect. 3. Implications of the revised $`T_{\mathrm{eff}}`$ scale and remaining uncertainties are discussed in Sect. 4.

## 2 Model ingredients

We have constructed spherically expanding non-LTE line-blanketed model atmospheres using the CMFGEN comoving-frame code of Hillier & Miller (hm98 (1998)). This code solves the equations of statistical equilibrium, radiative transfer, and radiative equilibrium, and allows for a direct treatment of line blanketing through the use of a super-level approach. The following ions are included in our calculations: H, He i-ii, C ii-iv, N ii-v, O ii-vi, Si ii-iv, S iv-vi, and Fe iii-vii, whose $``$ 2000 levels are described by $``$ 700 super-levels, corresponding to a total of $``$ 20000 bound-bound transitions.

For simplicity a constant Doppler profile (thermal width corresponding to the mass of Helium and $`T=20000`$ K plus a microturbulent velocity of $`v_{\mathrm{turb}}=20`$ km/s) is assumed for all lines in the statistical equilibrium and radiative transfer computation. To examine if a constant thermal width and the use of the large microturbulent velocity does not artificially enhance the photospheric blanketing, we have made test calculations with the correct depth and ion dependent thermal width and $`v_{\mathrm{turb}}=0.1`$ km/s. No significant changes in atmospheric structure, level populations, and the emergent spectrum were found. This is explained in part by the high density of lines in the UV part of the spectrum, which implies an average spacing between lines which is smaller than the typical Doppler width. The opacity in the wing of a line is therefore mostly dominated by the core opacity of the neighbouring line, and the exact intrinsic line profile is of little importance. With our standard choice, $``$ 80000 frequency points are necessary to correctly sample all lines.

The input atmospheric structure, connecting smoothly the spherically extended hydrostatic layers with the wind (parametrised by the usual $`\beta `$-law), is calculated as in Schaerer & de Koter (sdk (1997)) with the ISA-WIND code of de Koter et al. (dk96 (1996)) As the approximate temperature structure in ISA-WIND differs from the final radiative equilibrium temperature structure, the atmosphere structure in the quasi-hydrostatic part may be inconsistent with the final gas pressure gradient. However, for the issues discussed here the differences are small (corresponding to a change of $``$ 0.1 dex in $`\mathrm{log}g`$). In any case, the lines considered here are formed in the transition region whose structure/dynamics remain largely parametrised. The formal solution of the radiative transfer equation yielding the detailed emergent spectrum allows for incoherent electron scattering and includes standard Stark broadening tables for H, He i, and He ii lines. Our standard calculations assume $`v_{\mathrm{turb}}=5`$ km/s.

We have computed a grid of models representative of O dwarfs in the temperature range between $``$ 30000 and 50000 K. The model parameters are taken from the CoStar models A2-E2 of Schaerer & de Koter (sdk (1997)), with an additional model Y2 at ($`T_{\mathrm{eff}}`$,$`\mathrm{log}g`$) $``$ (31500,4.0) and the remaining parameters 111$`M=16.83`$ M, $`\mathrm{log}T_{\mathrm{eff}}=4.498`$, $`\mathrm{log}(L/L_{})=4.552`$, $`R=6.358`$ R, $`\mathrm{log}\dot{M}=7.204`$ M/yr, and $`v_{\mathrm{}}=2500`$ km/s. taken from stellar tracks of Meynet et al. (meynet94 (1994)). For each parameter set a line blanketed model with solar metallicity and a pure H-He model was computed.

## 3 Results

### 3.1 Blanketing effect on the temperature scale

The optical He i $`\lambda `$4471 and He ii $`\lambda `$4542 classification lines are used to assign spectral types to our models. Fig. 1 shows the effective temperature as a function of $`\mathrm{log}W^{}\mathrm{log}W(4471)\mathrm{log}W(4542)`$ and the corresponding spectral type according to Mathys (mat88 (1988)). The pure H-He models (open circles) follow closely the $`T_{\mathrm{eff}}`$–scale for dwarfs of VGS96, which is based on a compilation of stellar parameters determined using pure H-He plane parallel non-LTE model atmospheres. The comparison shows that if we neglect line blanketing our dwarf model grid would yield nearly the same absolute $`T_{\mathrm{eff}}`$–scale as the pure H-He plane parallel models adopted for the spectral analysis included in the compilation of VGS96.

The line blanketed model sequence (Fig. 1, filled symbols) shows a systematic shift to earlier spectral types for a given temperature, or equivalently a shift to lower $`T_{\mathrm{eff}}`$ for line blanketed models at a given spectral type. The difference ranges from $``$ 1500 K at spectral type O9.5 to $``$4000 K at spectral type O3 (cf. Fig. 1, solid line in lower panel). The difference with the VGS96 scale is shown as the dotted line. Our line blanketed scale smoothly joins earlier calibrations at O9.7V (see VGS96, Fig. 1).

As a spectral type corresponds to a given ionisation state of Helium in the line formation region, blanketed models must be more ionised than unblanketed models. The introduction of line blanketing leads to three main effects illustrated in Fig. 2 for the case of model C2 (cf. Figs. 13 and 14 of Schaerer & Schmutz ss94 (1994)). Qualitatively the same trends are obtained for all models.

* Blanketing leads to the backscattering of photons towards the inner atmosphere which forces the local temperature to rise so that flux conservation is fulfilled (backwarming effect; see upper panel).
* At the same time the radiation field becomes more diffuse, as quantified by the dilution factor $`\stackrel{~}{W}=1\frac{1}{4}F/J`$ shown in the middle panel, causing an increase of the mean intensity (cf. Abbott & Hummer ah85 (1985), Schaerer & Schmutz ss94 (1994)).
* In the outer part of the atmosphere ($`\mathrm{log}\tau _{\mathrm{Ross}}`$ –2 in Fig. 2) the ionisation is essentially controlled by the EUV flux, which is quite strongly reduced due to the blocking by numerous metal lines shown in Fig. 3. Here this effect dominates over 2), in contrast with the finding of Schaerer & Schmutz (ss94 (1994)), leading to a lower ionisation.

Effects 1) and 2) lead to a higher ionisation in the formation region of the classification lines. This results predominantly in an increase of W(4542) at $`T_{\mathrm{eff}}`$ $``$ 38000 K and a decrease of W(4471) at higher $`T_{\mathrm{eff}}`$(cf. Fig. 4).

Given the stronger mass loss and the corresponding increase of the wind density, one expects even larger temperature differences between non-blanketed and line blanketed models for giant and supergiant luminosity classes (cf. Abbott & Hummer ah85 (1985), Schmutz s98 (1998), Crowther et al. cro01 (2001)).

### 3.2 Dependence on model and stellar parameters

How strongly do our results depend on poorly known parameters such as the velocity law in the photosphere-wind transition zone, $`v_{\mathrm{turb}}`$, and variations of gravity and $`\dot{M}`$ expected within the dwarf class ? Do our calculations still miss opacity sources ?

As pointed out by Schaerer & Schmutz (ss94 (1994)) changes in He line profiles due to modifications of the velocity law $`v(r)`$ in the photosphere–wind transition zone can lead to similar equivalent widths variations as line blanketing. Test calculations for models A2 and C2 varying the slope $`\beta `$ from 0.8 (our standard value) to 1.5 show that both H-He and line blanketed models exhibit a similar shift in $`\mathrm{log}(W^{})`$. The obtained relative $`T_{\mathrm{eff}}`$ difference between H-He and blanketed models remains thus identical. The blanketed models with $`\beta =1.5`$ have $`\mathrm{log}(W^{})`$ lowered by $``$ 0.1–0.2 dex. However, as H$`\alpha `$ profile fits for O dwarfs are generally quite compatible with $`\beta 0.8`$ (e.g. Puls et al. puls96 (1996)), we do not expect drastic changes of the absolute scale from this effect.

An increase of the microturbulent velocity $`v_{\mathrm{turb}}`$ from 5 to 20 km/s in blanketed models increases the strength of He i $`\lambda `$4471 (cf. Smith & Howarth sh98 (1998), Villamariz & Herrero vh00 (2000)), and leads to a shift of $``$ \+ 0.05 to 0.1 dex in $`\mathrm{log}(W^{})`$ (i.e. towards later types) for models with $`T_{\mathrm{eff}}`$ $``$ 42000 K. For hotter stars the difference is negligible.

The effect of line blanketing is strengthened further in denser winds (cf. Abbott & Hummer ah85 (1985), Schmutz s98 (1998)). Models C2 and D2 with an increased mass loss rate by a factor of 2 show a shift of $`\mathrm{log}(W^{})`$ between $``$ – 0.05 and – 0.1 dex.

Test calculations for model C2 including also Nickel (Ni iv-vi) show unchanged He lines. Other models including also Ar, Ne, and Ca confirm that Fe blanketing dominates.

While microturbulence and mass loss affect (though in opposite ways) the exact $`T_{\mathrm{eff}}`$-scale, their exact importance will have to be studied in future comparisons.

### 3.3 Comparison with observations

As a first comparison of our models with observations we show in Fig. 4 the predicted and observed equivalent widths of He i and He ii classification lines and other strong He lines frequently used in spectral analysis. The observational data is taken from Mathys (mat88 (1988), mat89 (1989)) and Conti & Alschuler (ca71 (1971)). The observational scatter is real, as the typical measurement errors are $``$ 5–7 %. The general trend is that the He i $`\lambda `$4471 and He i $`\lambda `$4388 equivalent widths are well represented by the models, while He ii $`\lambda `$4542 seems to be overestimated by $``$ 20% for spectral types earlier than O7. He ii $`\lambda `$4200 behaves as He ii $`\lambda `$4542 . The other equivalent widths remain essentially unchanged by all other parameter variations discussed above (Sect. 3.2). A value of $`\beta 1.5`$, a stronger increase of $`\dot{M}`$, or an unrealisticly large reduction of $`\mathrm{log}g`$ would be necessary to reduce the predicted equivalent widths of the Stark broadened He ii lines.

Strictly speaking, if we were simply to reduce $`W(4542)`$ by $``$ 20% while keeping $`W(4471)`$ constant for early spectral types, this would result in a change of $``$ – 0.08 dex in $`\mathrm{log}(W^{})`$ thus reducing the shift in the $`T_{\mathrm{eff}}`$–scale between line blanketed and pure H-He models from $``$ 4000 K to 3000 K in the high temperature part. Future tailored spectral analysis should allow to assess more precisely the achievable fit accuracy and the precise importance of the parameters discussed in Sect. 3.2 on the stellar parameters.

### 3.4 Comparison with previous analysis

As discussed in Sect. 1, few earlier studies have addressed the effect of line blanketing in O stars. Essentially all investigations concur with a reduction of $`T_{\mathrm{eff}}`$ when blanketing is included.

Abbott & Hummer (ah85 (1985)) have constructed a core-halo model where backscattered radiation due to multiple line scattering in the wind modifies the plane parallel photosphere. Their so-called “wind blanketed” models yield a decrease of $`T_{\mathrm{eff}}`$ by $``$ 10 % for O4 types (similar to our results), $``$ – 2000 K for an O9.5 supergiant, but essentially no shift for O9.5 dwarfs (Bohannan et al. bo90 (1990), Voels et al. vo89 (1989)). The latter finding is likely due to lack of photospheric blanketing (inherent to their method) and modest wind blanketing due to the comparatively low mass loss rates of O9.5 dwarfs.

An improved Monte-Carlo opacity sampling method of a unified photosphere–wind model was used by Schaerer & Schmutz (ss94 (1994)), Schaerer & de Koter (sdk (1997)), and subsequently applied to a larger parameter space by Schmutz (s98 (1998)). For mass loss rates comparable to the values adopted here (typical for dwarfs with low mass loss) the models of Schmutz (s98 (1998)) indicate differences from $``$ – 600 K at O8 to $``$ – 2000 K at O4, which is half the shift deduced from Fig. 1 and roughly the difference obtained with $`Z=1/8`$ Z (see Sect. 4). This indicates that their method underestimates line blanketing compared to CMFGEN.

Using plane parallel line blanketed non-LTE models based on opacity distribution functions Hubeny, Heap & Lanz (hub98 (1998)) found that a pure H-He model with $`T_{\mathrm{eff}}`$ $``$ 37500 K and $`\mathrm{log}g`$=4.0 is necessary to reproduce the H and He lines of a line blanketed model with $`T_{\mathrm{eff}}`$ = 35000 K and same gravity. As can be seen from Fig. 1 our results are in excellent agreement with their result.

LTE line blocking has been included in plane parallel models by Herrero et al. (her00 (2000)) primarily to resolve discrepancies between He i singlet and triplet lines. For stars with $`T_{\mathrm{eff}}`$ $``$ 40000 K this leads to a strengthening of He i $`\lambda `$4471, opposite to the effect found in all above studies including ours. This results must be due to an incomplete treatment of the various effects of line blanketing (cf. above), and appears to be unphysical. This discrepancy with line blanketed models has also been noted by the authors.

## 4 Implications and concluding remarks

The importance of line blanketing obviously depends on metallicity $`Z`$. Therefore one may wonder at which $`Z`$ the stellar parameters will again correspond to the results obtained with pure H-He (metal-free) atmosphere models, i.e. close to the VGS96 scale. Test calculations for models A2 and D2 with a metallicity close to the SMC value (1/8 Z) show still a reduction of $`T_{\mathrm{eff}}`$ compared to pure H-He models : $`\mathrm{\Delta }`$ $`T_{\mathrm{eff}}`$ is $``$ 60 % that found at solar metallicity.

As the bolometric correction is essentially unchanged by line blanketing, and the $`M_V`$ versus spectral type ($`Sp`$) calibration independent to first order from modeling, we can use the BC-$`T_{\mathrm{eff}}`$ relation of VGS96 to derive luminosities through $`\mathrm{log}(L/L_{})=2.736\mathrm{log}T_{\mathrm{eff}}(Sp)0.4M_V(Sp)9.164`$. This relation shows that the predicted reduction of $`T_{\mathrm{eff}}`$ by $``$ 0.04 dex implies a downward revision of $`L`$ by $``$ 0.1 dex for dwarfs of a given spectral type.

Since line blanketing is mostly efficient in the EUV, the ionising spectrum below 912 Å is modified. The total number of Lyman continuum photons $`Q_0`$ predicted by our models is in good agreement with the calculations of Schaerer & de Koter (sdk (1997)). The change of $`Q_0`$ due to the shift in the $`T_{\mathrm{eff}}`$-Sp calibration, taking into account the change of both the radius and the ionising flux per unit surface area $`q_0`$, is given by $`\mathrm{\Delta }(\mathrm{log}Q_0)=1.264\mathrm{\Delta }(\mathrm{log}T_{\mathrm{eff}})+\mathrm{\Delta }(\mathrm{log}q_0(T_{\mathrm{eff}}))`$, where the latter term is dominant (see Schaerer & de Koter sdk (1997)). For a given spectral type between O4V and O9V this amounts typically to a reduction of $`Q_0`$ by $``$ 40 %.

While the results presented here provide a clear improvement over earlier calibrations, and a general reduction of $`T_{\mathrm{eff}}`$ due to line blanketing is unavoidable, we wish to caution that the absolute $`T_{\mathrm{eff}}`$ scale may still be subject to revisions for the following reasons. First, tailored multi-wavelength analysis of individual objects are required to test the present models in more depth for O stars, as recently started by Bouret et al. (bouret\_smc (2001)), Hillier et al. (john (2001)), and Crowther et al. (cro01 (2001)). Second, the effect of X-rays on the overall ionisation balance and in particular on the Helium lines remains to be studied. Indeed for late O and B stars, depending on the relative X-ray to photospheric flux at energies close to the relevant ionisation potentials and the wind density, X-ray emission (likely due to shocks) is expected to increase the ionisation of most ions (MacFarlane et al. mac94 (1994)). Nonetheless, first test calculations with CMFGEN seem to indicate that photospheric lines are not affected by X-rays generated in the wind. Finally, we note that comparisons of photoionisation models calculated using fluxes from recent atmosphere models (including CMFGEN and Pauldrach et al. p01 (2001) models) with ISO observations of H ii regions possibly reveal a flux deficiency at energies $``$ 34.8–40.9 eV (Morisset et al. mor01 (2001), but cf. Giveon et al.giv01 (2001)). The importance of the latter two findings — possibly related to each other — on the lines used here as $`T_{\mathrm{eff}}`$ indicators remains to be studied.

As UV and optical classification lines of O stars depend in fact on several parameters ($`T_{\mathrm{eff}}`$, gravity, mass loss rate, metallicity, rotation; e.g. Abbott & Hummer ah85 (1985), Schmutz s98 (1998), Walborn et al. wal95 (1995)), spectral type and luminosity class calibrations must ultimately account for this multi dimensionality. Some of these issues will be addressed in subsequent publications.

###### Acknowledgements.

We thank the “Programme National de Physique Stellaire” (PNPS) for support for this project and the CALMIP center in Toulouse for generous allocation of computing time.


text:
# Statistical mechanics of screened spatially indirect excitons

## I Introduction

Spatially separated quasi-two-dimensional (2D) electron-hole plasma can be found in different modern nano-heterostructure-based devices. There is a continuous attention to closely situated (so called double) quantum wells. This attention is mostly due to the possibility of the Bose-Einstein condensation of spatially indirect excitons in such a system. Another physical object with spatially separated electron-hole plasma, which currently attracts increasing interest, is a GaN/AlGaN quantum well, grown on a sapphire substrate. In the latter case electrons and holes are separated by huge (of the order of one MV/cm) internal electric field. Because of the enhanced binding energy in wide-gap semiconductors, excitons in nitride-based nanostructures are a subject of extensive experimental and theoretical research.

The aim of this work is to provide a consistent many-body theory, describing statistical mechanics of the quasi-2D exciton/free carrier plasma with spatial charge separation.

## II Screening in the Spatially-Separated Electron-Hole System

First of all, the potentials of electron-electron and electron-hole interaction should be established. The Thomas-Fermi (or, as it is often called in the 2D case, Stern-Howard) approximation is one of the most frequently used methods for dealing with static screening of Coulomb interaction in semiconductor systems. In the frame of this approximation the charge, induced by external perturbation, originates from a redistribution of free carriers and is proportional to a local electric field. The resulting electric field potential is given by the Poisson equation in the form:

$$\mathrm{\Delta }\varphi 2q_s^e\overline{\varphi _e}g_e(z)2q_s^h\overline{\varphi _h}g_h(z)=\frac{4\pi }{\kappa }Q_{ext},$$
(1)
where $`\overline{\varphi _a}(\rho )=\varphi (\rho ,z)g_a(z)𝑑z`$, $`g_a(z)`$ is a free-carrier-distribution function in the direction normal to the quantum well plane, $`\kappa `$ is a background dielectric constant, $`Q_{ext}`$ is an external charge density, and the screening parameter in the 2D case is

$$q_s^a=\frac{2\pi e^2}{\kappa }\frac{n_a}{\mu _a}=\frac{2}{a_B}\frac{m_a}{m_{eh}}\left[1\mathrm{exp}\left(\beta \pi \mathrm{}^2n_a/m_a\right)\right].$$
(2)
Here $`n_a`$ is a free particle density, $`\mu _a`$ is a quasi-Fermi level, $`m_a`$ is an effective mass, $`\mu _{eh}=m_em_h/(m_e+m_h)`$ is the reduced effective mass, $`a_B=\kappa \mathrm{}^2/(\mu ^{}e^2)`$ is the 3D exciton Bohr radius, and $`\beta =1/(k_BT)`$.

The contribution of excitons to screening in the conventional single quantum well can be treated within the Thomas-Fermi approximation as a reduction of free carrier densities entering Eq.(2). In the case of spatially separated plasma, screening by excitons can not be neglected even in the Thomas-Fermi approximation. The indirect exciton has the finite dipole moment, and the external electric field causes redistribution of excitons, which plays its own role in the resulting field reduction. Using the same approach as for a free-carrier plasma, the indirect-exciton part of the induced charge density can be written as

$$Q_{ind}^{exc}=\frac{\kappa }{2\pi }q_s^{exc}(g_h(z)g_e(z))\left[e\overline{\varphi _h}(\rho )e\overline{\varphi _e}(\rho )\right].$$
(3)
Here we introduced the excitonic screening parameter

$$q_s^{exc}=\frac{2\pi e^2}{\kappa }\frac{n_{exc}}{\mu },$$
(4)
where $`n_{exc}`$ is the exciton density and $`\mu =\mu _e+\mu _h`$.

In the present paper we consider a simple model of spatially separated electron-hole plasma, where electrons and holes are purely two-dimensional and reside in different planes with distance $`d`$ between them. The generalization of this model to incorporate finite-width quantum wells will be presented somewhere else. The two-plane model is particularly valuable for studying the interplay between screening and spatial separation of different types of carriers. The problem of finding electron-electron repulsion and electron-hole attraction potentials has cylindrical symmetry, so we use the following representation:

$$\varphi (\rho ,z)=\frac{1}{2\pi }_0^{\mathrm{}}\varphi _q(z)J_0(\rho q)q𝑑q,$$
(5)
where $`J_0`$ is the Bessel function. Let us assume that the electronic charge resides in the plane with coordinate $`z=d`$ and the holes are in the plane $`z=0`$. In this model the distribution functions $`g_a(z)`$ are given by delta-functions. The Poisson equation is reduced to a one-dimensional problem:

$$\frac{d^2}{dz^2}\varphi _q(z)q^2\varphi _q(z)2q_s^h\varphi _q(0)\delta (z)2q_s^e\varphi _q(d)\delta (zd)$$

$$2q_s^{exc}(\delta (z)\delta (zd))[\varphi _q(0)\varphi _q(d)]=\frac{4\pi e}{\kappa }\delta (zd)$$
(6)
Solving Eq.(6), we obtain repulsion and attraction potentials in $`q`$-representation:

$$V_s^{eh}(q)=e\varphi _q(0)=\frac{2\pi e^2}{\kappa }\frac{qe^{qd}+q_s^{exc}(1e^{2qd})}{\mathrm{\Theta }+q_s^{exc}\mathrm{\Xi }},$$
(7)
$$V_s^{ee}(q)=e\varphi _q(d)=\frac{2\pi e^2}{\kappa }\frac{q+(q_s^h+q_s^{exc})(1e^{2qd})}{\mathrm{\Theta }+q_s^{exc}\mathrm{\Xi }},$$
(8)
where both terms $`\mathrm{\Xi }=[1\mathrm{exp}(qd)][2k+(q_s^e+q_s^h)(1+\mathrm{exp}(qd))],`$ and $`\mathrm{\Theta }=(q+q_s^{(e)})(q+q_s^{(h)})q_s^{(e)}q_s^{(h)}\mathrm{exp}(2qd)`$ in denominators of Eqs.(7,8) do not depend on the excitonic screening parameter $`q_s^{exc}`$.

## III Correlated density

In order to provide a consistent treatment of many-body effects in a 2D electron-hole plasma, we base our theory on a quasi-particle picture and Green functions technique. Electrons and holes are redistributed among the quasi-particle states with a sharp spectral peak and dispersion $`ϵ_a=\mathrm{}^2k^2/2m_a`$ (where the energy is measured from the corresponding size-quantized level renormalised by the many-body shift) and some other ”correlated” states originating from the interaction between quasi-particles.

The total density of electrons can be divided into two parts, the first part representing “free” particles with a quasi-Fermi level $`\mu _e`$, and the second part incorporating interaction between quasiparticles. This second part, which is called the correlated density, is given by

$$n_e^{corr}=2\underset{k}{}\frac{d\mathrm{}\omega }{\pi }\mathrm{\Gamma }_e(k,\omega )\left[f_e(\mathrm{}\omega )f_e(ϵ_e(k))\right]\frac{}{\mathrm{}\omega }\frac{\mathrm{P}}{ϵ_e(k)\mathrm{}\omega }.$$
(9)
Here $`\mathrm{\Gamma }_e(k,\omega )=\mathrm{Im}\mathrm{\Sigma }_e(k,\omega )`$ is an imaginary part of the single-particle self-energy, $`f_a(ϵ)=[\mathrm{exp}(\beta (ϵ\mu _a))+1]^1`$ is the Fermi distribution function, and the derivative of the principal value is given

$$\frac{}{\mathrm{}\omega }\frac{\mathrm{P}}{ϵ\mathrm{}\omega }=\underset{\nu 0}{lim}\frac{(ϵ\mathrm{}\omega )^2\nu ^2}{[(ϵ\mathrm{}\omega )^2+\nu ^2]^2}.$$

From the point of view of diagrammatic technique, exciton is a succession of electron-hole scatterings, which is represented by the series of diagrams serving as a basis for the ladder approximation. We use the ladder approximation to our many-body problem, keeping only ladder-type diagrams in the expression for the electron self energy:

$$\mathrm{\Sigma }_e(\stackrel{}{k},\mathrm{\Omega }_e)=\frac{1}{\beta }\underset{b,\stackrel{}{k}^{},z_b}{}\left[2\stackrel{~}{T}_{eb}(\stackrel{}{k},\stackrel{}{k}\stackrel{}{k}^{},\stackrel{}{k},\mathrm{\Omega }_e+\mathrm{\Omega }_b)\delta _{eb}\stackrel{~}{T}_{eb}(\stackrel{}{k}^{},\stackrel{}{k}\stackrel{}{k}^{},\stackrel{}{k},\mathrm{\Omega }_e+\mathrm{\Omega }_b)\right]G_b(\stackrel{}{k},\mathrm{\Omega }_b),$$
(10)
where the summation is taken over different types of particles, $`b=e,h`$, 2D wave-vector $`\stackrel{}{k}^{}`$ and Fermi-type Matsubara frequencies $`\mathrm{\Omega }_b`$. The quantity $`\stackrel{~}{T}`$ represents a sum of ladder diagrams for electron-electron or electron-hole scattering and is given by the frequency-dependent T-matrix equation:

$$\stackrel{~}{T}_{ab}(\stackrel{}{k},\stackrel{}{q},\stackrel{}{k}^{\prime \prime },\mathrm{\Omega })=$$

$$V_s^{ab}(\stackrel{}{k}\stackrel{}{k}^{\prime \prime })+\underset{\stackrel{}{k}^{}}{}V_s^{ab}(\stackrel{}{k}\stackrel{}{k}^{})G_{ab}(\stackrel{}{k}^{},\stackrel{}{q},\mathrm{\Omega })\stackrel{~}{T}_{ab}(\stackrel{}{k}^{},\stackrel{}{q},\stackrel{}{k}^{\prime \prime },\mathrm{\Omega }),$$
(11)
where the potentials $`V_s^{ab}`$ are given by Eqs.(7,8),

$$G_{ab}(\stackrel{}{k},\stackrel{}{q},\mathrm{\Omega })=\frac{N_{ab}(\stackrel{}{k},\stackrel{}{q})}{ϵ_a(\stackrel{}{k})+ϵ_b(\stackrel{}{q}\stackrel{}{k})\mathrm{}\mathrm{\Omega }},$$

and

$$N_{ab}(\stackrel{}{k},\stackrel{}{q})=1f_a(ϵ_a(\stackrel{}{k}))f_b(ϵ_b(\stackrel{}{q}\stackrel{}{k})).$$

Following Zimmermann , who derived a similar formula for the 3D case, we obtain the expression for the correlated density in the 2D electron-hole system:

$$n_a^{corr}=\frac{2}{\beta \pi \mathrm{}^2}\underset{m=\mathrm{}}{\overset{+\mathrm{}}{}}\left[\underset{n}{}M_{eh}L_{eh}\left(ϵ_{mn}\right)+\frac{1}{\pi }\underset{b=e,h}{}M_{ab}\lambda _m^{ab}_0^{\mathrm{}}𝑑kL_{ab}\left(\frac{\mathrm{}^2k^2}{2m_{ab}}\right)2\mathrm{sin}^2\delta _m^{ab}(k)\frac{d\delta _m^{ab}}{dk}\right],$$
(12)
$$L_{ab}(ϵ)=\mathrm{ln}\left[1\mathrm{exp}\left(\beta (\mu _a+\mu _bϵ)\right)\right],$$

$$M_{ab}=m_a+m_b,m_{ab}=\frac{m_am_b}{m_a+m_b},\lambda _m^{ab}=1\delta _{ab}(1)^m/2.$$

Here $`ϵ_{mn}`$ are the values of the exciton bound state energies, where the first index $`m`$ denotes the angular momentum and the second index $`n`$ numerates different states for a given $`m`$. The bound state energies correspond to the poles of the T-matrix. The quantities $`\delta _m^{ab}(k)`$ give the phases of the coefficients in the Fourier series expansion of the T-matrix. The values of these quantities coincide with the scattering phase shifts when the k-space filling is neglected. Thus, the second part of Eq.(12) describes the scattering states contribution to the correlated density. In the low density/high temperature limit this expression differs from that of Portnoi and Galbright  only by the factor $`2\mathrm{sin}^2\delta _m^{ab}`$ in the integrand. The “e-e” part of the correlated density has a negative sign and is responsible for the reduction of the plasma density due to repulsion between the carriers of the same type. The electron-hole part of the density is positive and can be attributed to the presence of the excitons. As it was shown in Ref. for dilute plasma, this part of the density does not change abruptly when one of the bound states disappears with increasing screening. We treat the electron-hole part of the correlated density as the density of excitons:

$$n_{exc}=\frac{2}{\beta }\frac{M_{eh}}{\pi \mathrm{}^2}\underset{m=\mathrm{}}{\overset{+\mathrm{}}{}}\left[\underset{n}{}L_{eh}\left(ϵ_{mn}\right)+\frac{1}{\pi }_0^{\mathrm{}}𝑑kL_{eh}\left(\frac{\mathrm{}^2k^2}{2m_{eh}}\right)2\mathrm{sin}^2\delta _m^{eh}(k)\frac{d\delta _m^{eh}}{dk}\right].$$
(13)
Substituting $`n_{exc}`$ from Eq.(13) into Eq.(4) we get the following expression for the excitonic screening parameter:

$$q_{exc}=\frac{4M_{eh}e^2}{\kappa \mathrm{}^2}\underset{m=\mathrm{}}{\overset{+\mathrm{}}{}}\left[\underset{n}{}f_B\left(ϵ_{mn}\right)+\frac{1}{\pi }_0^{\mathrm{}}𝑑kf_B\left(\frac{\mathrm{}^2k^2}{2m_{eh}}\right)2\mathrm{sin}^2\delta _m^{eh}(k)\frac{d\delta _m^{eh}}{dk}\right],$$
(14)
$$f_B(ϵ)=\left[\mathrm{exp}(\beta (ϵ\mu _e\mu _h))1\right]^1.$$

To calculate the exciton density we have to run the following self-consistent procedure. Firstly, for given electron-hole density we calculate the screening parameters from Eq.(2) and obtain corresponding bound state energies and phase shifts. Then, using these results we calculate the new values of excitonic screening parameter and free carrier density and recalculate bound state energies and phase shifts, repeating this procedure until the exciton density stops changing.

## IV Results and discussion

Throughout the calculations we neglect the phase shift dependence on the carrier density and obtain the values of $`\delta _{ab}`$ using the 2D modification of the variable phase approach. To find the bound state energies we use the method, proposed by Campi and co-workers, which is also based on the ladder approximation. We take into account only the ground exciton state. The material parameters are chosen to match these of the GaAs/AlGaAs quantum wells, and the difference between dielectric constants of the QW and the barrier is neglected.

In Figure 1 we present the dependence of the binding energy of the exciton ground state on the inter-plane distance. One can see (Fig. 1b) that the decrease of the exciton binding energy with increasing the distance is substantially enhanced by the presence of free carriers and indirect excitons, which screen the Coulomb attraction. Figure 2 shows the exciton density as a function of the total plasma density for different inter-plane distances. One can see that the exciton density has a maximum and then decreases sharply with increasing the total number of carriers, which reflects the action of both screening and k-space filing.


text:
# Work functions of pristine and alkali-metal intercalated carbon nanotubes and bundles

## Abstract

The work functions (WF) of single-walled carbon nanotubes and bundles are studied using first principles methods. For individual metallic tubes, the WF is independent of the chirality and increase slightly with tube diameter. For semiconducting tubes, the WF (as defined by the HOMO energy) decreases rapidly. The WF of nanotube bundles ($``$ 5 eV) shows no clear dependence on the tube size and chirality, slightly higher than individual tubes. Calculations on finite tubes show no substantial difference in the tube end and the side wall. Upon alkali-metal intercalation, the WF decreases dramatically and the electronic states near the Fermi level are significantly modified. The metallic and semiconducting nanotubes bundles become indistinguishable.

The work function is one of the critical quantities in understanding the field emission properties of carbon nanotubes . Although the work function can be estimated from the field-emission spectra based on Fowler-Nordheim model , the results are not reliable due to the uncertainty of the local geometry of nanotubes . Other experimental measurements on work functions of both single-walled nanotubes (SWNTs) and multi-walled nanotubes (MWNTs) included ultraviolet photoemission spectroscopy (UPS)  and transmission electron microscopy . It was found that the work functions of MWNTs are about 0.1-0.2 eV lower than that of the graphite , while the SWNT bundles have slightly higher work functions . Upon Cs intercalation the work functions of carbon nanotubes are reduced dramatically , which leads to a significant enhancement in field emission . Up to now, theoretical works were only limited to a few nanotubes with finite-lengths . In this letter, we report results of first principles calculations on the work functions of individual carbon nanotubes and nanotube bundles. The effect of tube diameter and chirality, and the difference between the tube end and the side wall are investigated. The alkali-metal (K, Rb, Cs) intercalated carbon nanotube bundles are also discussed.

The work function of a bulk metal is related to its Fermi energy $`E_F`$ by $`WF=\varphi E_F`$, where $`\varphi `$ is the electrostatic potential caused by “spilling out” of electron density at the metal surface . For those metals with low electron density such as K, Rb, Cs, it is known that $`\varphi `$ is much smaller than $`E_F`$ . For carbon nanotubes, the conduction electron density is much smaller than that of K or Cs . In this work, we approximate the WF by the Fermi energy $`E_F`$ . To determine the Fermi level of carbon nanotubes with respect to the vacuum level, we perform all electron LCAO calculations based on the density functional DMol package . The density functional is treated by the generalized gradient approximation (GGA)  with the exchange-correlation potential parameterized by Wang and Perdew (PW91) . For the infinite nanotubes, one-dimensional (1-D) periodic boundary condition is applied along the tube axis. The Brillouin zone is sampled by large sets of Monkhorst-Pack k meshes  (Along the tube axis, 40 k points are used for standard calculations and 160 points for accurate electronic density of states). Benchmark calculations are carried out on several alkali metals solids and the graphene sheet. The ionization potential of a C60 cluster is also calculated from the total energy difference between the neutral and the charged cluster. As shown in Table I, the overall agreement between theory and experiments  for different systems is reasonable.

We first address the work functions of the individual metallic SWNTs of infinite length, with diameter ranging from 5.7Å to 16.3 Å . Both the armchair ($`n`$,$`n`$) and the zigzag ($`n`$,0) ($`n=3m`$) chirality are considered. The calculated WFs range from 4.63 to 4.77 eV, given in Fig.1 and Table II. Careful examinations show that the WF of the metallic tubes increases slightly with 1/D (Fig.1). Extrapolation towards the D$`\mathrm{}`$ limit gives a WF of 4.83 eV, close to the calculated value for the graphene sheet (4.8 eV). It is worthy to note that the WFs of both the armchair and the zigzag metallic tubes fit the same linear dependence, indicating that the WFs of metallic nanotubes are independent of chirality.

For the purpose of comparison, we define the WF of semiconducting tube as the highest occupied molecular orbital (HOMO) energy. The WFs of the semiconducting tube are substantially higher than the metallic ones (Fig.1). The calculated WF decreases linearly with 1/D and approaches an extrapolation limit 4.73 eV at D$`\mathrm{}`$. The strong diameter dependence come from the well-known gap decrease with tube diameter in semiconducting tubes . There is no direct experiment measurement on the work functions of individual SWNTs. But recent experiments on the MWNT tips suggest that the WFs of semiconducting tubes ($`5.6`$ eV) are high than those of metallic ones ($`4.64.8`$ eV) .

To explore the difference in work function between the side wall and the tube end, we carry out model calculations on several finite (5,5) tubes with length up to 21.8 Å. As shown in Fig.2, one end of the nanotube is capped by a half-C60 while the other end is H-terminated. The calculated HOMO energy decreases with the tube length $`l`$ and approaches the infinite limit rapidly (4.87 eV, 4.63 eV, 4.55 eV for $`l=3.8\AA ,11.7\AA ,21.8\AA `$ respectively). In Fig.3, we plot the HOMO wave function distributions for the longest tube studied (160 carbon atoms on side wall plus 30 carbon atom on the end). The HOMO are clearly delocalized over the whole tube. The wave function at the tube end is comparable to that on the side wall. Similar distributions can be found for the molecular orbitals just below the HOMO level. Therefore, we expect that the work function at the tube end should not be much different from that of the sidewall. For the finite size nanotubes investigated, we find interesting oscillations in both the HOMO energy and the wave function distribution with the tube length. The oscillation period is three times that of the unit cells length along the tube axis, corresponding to the Fermi wavelength of armchair tubes . This oscillation behavior implies that care should be taken when extrapolating the finite size results to the infinite limit . Our detailed investigations on the electronic properties of finite nanotubes with different lengths and ends will be presented elsewhere.

We now discuss the work functions of the nanotube bundles. The bundles are modeled by two dimensional hexagonal lattices of uniform SWNTs . For all the metallic tube bundles studied, the calculated work functions are around 5 eV (see Table II), slightly higher than individual metallic nanotubes and the graphene. The increase of WFs in tube bundles can be understood by the tube-tube interaction . From the experiments based on UPS, the work functions of SWNT bundles are determined to be 4.8 eV  or 5.05 eV . Our theoretical results agree well with the available data on nanotube bundles.

In addition to pristine materials, the electronic properties of carbon nanotubes can be efficiently controlled by alkali-metal intercalations. Recently, nanoelectronics devices have been constructed on the basis of alkali-metal doped carbon nanotubes . Thus, it is important to investigate the effect of alkali-metal intercalations. We have carried out systematical calculations on the effect of intercalations in both metallic (10,10) and semiconducting (17,0) tubes bundles. Intercalation density up to K0.1C (Rb, Cs) (close to the saturation density in nanotube bundles  and the graphite ) are studied. Both the inside of SWNTs and the interstitial sites are explored. The initial configurations of alkali-metal atoms are chosen to maximize the ion-ion distance . Structural relaxations find no significant change from the initial configurations. For a given intercalation concentration, the WF is insensitive to the detailed configurations of the intercalated atoms. The lattice constants of two-dimensional hexagonal lattices are expanded. For example, we find that the intercalations of K atoms into the interstial sites of (10,10) tube bundle can induce about $``$2Å lattice expansion, which is comparable to expansion of 1.95 Å in K intercalated graphite  and 1.85 Å in the HNO3 intercalated SWNT bundles . In Fig.3, we present the calculated work functions of intercalated tube bundles as functions of the metal concentration for various metals. In general, the work function dramatically decreases after small amount of alkali-metal intercalations. The reduction becomes much slower at higher intercalation density. Furthermore, there is almost no difference between the WFs of intercalated (10,10) and (17,0) tube bundles. Thus, the metallic and semiconducting tube bundles become indistinguishable. The WF for Cs-doped nanotube bundles is slightly lower than that of the Rb- or K- doped systems. This trend from K to Cs is consistent with the WFs for bulk metal that ($`WF(K)>WF(Rb)>WF(Cs)`$, see Table I).

Intuitively, the reduction of WFs can be understood by the charge transfer from metal to carbon nanotube, which shifts the EF of conduction band towards the vacuum level. Experimentally, the charge transfer from metal atoms to the nanotube have been confirmed by the resistance  and Raman spectrum . However, our present results show that interaction between K (or Rb, Cs) and nanotube cannot be simply described by a rigid band model with charge transfer. Shown in Fig.4 is the electronic density states for pristine and K-intercalated (10,10) SWNT bundles. Similar to the case of Li intercalation , the valence bands are almost not affected by K intercalations. In contrast, the conduction bands are significantly modified by potassium-carbon interaction. New peaks associated with alkali-metal atoms are found in the conduction bands. The density of states near the Fermi level is significant enhanced by the contributions of alkali metals. Our thoeretical results are supported by recent experiments on the optical properties of the K or Cs intercalated nanotube bundles . In addition, we have also calculated the density of states for intercalated semiconducting bundles. Similar to the case of work function, the metallic and semiconducting nanotubes bundles become indistinguishable. This is consistent with recent NMR measurements on K-intercalated SWNT bundles .

In summary, we have performed first principles calculations on the work functions of pristine and intercalated SWNT nanotube and bundles. The WFs of metallic nanotubes are weakly depend on the tube size and comparable to the graphene sheet. The work functions of all the metallic tube bundles are around 5 eV, slightly higher than those of individual tubes. Finite size tube calculations reveal the delocalized HOMO state, implying no substantial difference between the WF of the tube end and the side wall. Upon alkali-metal intercalation, the WFs of tube bundles decrease dramatically and the electronic states near Fermi level are significantly modified. The metallic and semiconducting nanotubes bundles become indistinguishable. This work is supported by the U.S. Army Research Office Grant No. DAAG55-98-1-0298, the Office of Naval Research Grant No. N00014-98-1-0597 and NASA Ames Research Center. The authors thank Prof. O.Zhou, Prof. Y.Wu for helpful discussions. We acknowledge computational support from the North Carolina Supercomputer Center.

E-mail address: zhaoj@physics.unc.edu
E-mail address: jpl@physics.unc.edu


text:
# A microscopic model for a class of mixed-spin quantum antiferromagnets

## Abstract

We propose a microscopic model that describes the magnetic behavior of the mixed-spin quantum systems R2BaNiO5 (R= magnetic rare earth). An evaluation of the properties of this model by Quantum Monte Carlo simulations shows remarkable good agreement with the experimental data and provides new insight into the physics of mixed-spin quantum magnets.

PACS numbers: 75.25.+z, 75.10.Jm, 75.50.Ee

J.V. Alvarez

Roser Valentí

A. Zheludev

Low-dimensional (low-D) quantum magnetism remains at the forefront of condensed matter research for almost two decades. This is primarily because this field deals with simple models of magnetism that demonstrate a broad spectrum of complex quantum-mechanical phenomena of a general and fundamental nature. The simplicity of the Hamiltonians involved, as well as the low dimensionality, often allow accurate theoretical and numerical treatment. On the other hand, the discovery of real low-D materials allows experimental studies for direct comparison with theory on the quantitative level. A particular area of great recent interest is the co-existence of “quantum” and “classical” properties in systems composed of weakly coupled quantum spin chains. One good example is the interplay between continuum and spin-wave dynamics in gapless quasi-one-dimensional (quasi-1D) Heisenberg systems such as KCuF3  and BaCu2Si2O7 . At the focus of the present work is a more complex phenomenon: the seemingly paradoxical co-existence of 3D magnetic long-range order and 1D quantum gap excitations in rare earth nickelates with the general formula R2BaNiO5 (R= magnetic rare earth).

R2BaNiO5 species have two types of spin carriers. $`S=1`$ Ni2+ ions form distinct antiferromagnetic chains running along the $`a`$ axis of the crystal structure. If these chains were perfectly isolated, they would be a classic example of a Haldane antiferromagnet with no long-range order even at $`T=0`$ and a gap in the magnetic excitation spectrum . Long-range ordering that occurs in R2BaNiO5 at low temperatures is driven by the magnetic $`s=\frac{1}{2}`$ R3+ ions, positioned in-between the chains. 3D magnetic order involves both the R3+ and Ni2+ spins. What is particularly interesting though, is that Haldane-gap excitations associated with the Ni-chains persist in the ordered state and co-exist with conventional order-parameter excitations (spin waves). By now, a wealth of experimental data on several R2BaNiO5 compounds has been accumulated . Most theoretical work, however, was based on a simple Mean Field (MF)  or Random Phase Approximation (RPA) treatment of the interactions between Ni-chain and rare earth subsystems. The existing numerical studies were aimed at calculating the “bare” properties of the Ni-chains, needed to complete the MF/RPA equations. While the MF/RPA approach turned out to be an extremely useful tool in understanding the basic physics involved, it has numerous intrinsic limitations, particularly for dilute (R1-xYx)2BaNiO5 systems , or in the case of strong Ni-R coupling. In the present work we abolish the MF/RPA framework and construct a microscopic model for R2BaNiO5 compounds. We then employ this model in a systematic first-principle numerical study of these interesting materials.

In order to construct the microscopic model we will use the available experimental results to design an appropriate description for the magnetic ions and interactions in the system: (i) A good understanding of the Ni-subsystem can be drawn from the known properties of Y2BaNiO5, a material where the magnetic rare earths have been replaced by non-magnetic Y3+. The Yttrium-nickelate is an almost perfect physical realization of a Haldane gap antiferromagnet. Neutron scattering data show a large gap $`\mathrm{\Delta }10`$meV and well characterized Haldane excitations. This suggests that the intrachain Ni-Ni exchange coupling in the chains running along the $`a`$-axis has to be $`J25`$meV. Furthermore, there is no evidence of 3D long range order even at the lowest temperatures. Any residual direct coupling between the chains, both in $`b`$ or $`c`$ direction, is smaller than the critical value necessary to establish AF order at low enough temperatures . This critical value depends quadratically on the gap and it is expected to be large in this family of nickelates due to the magnitude of the gap. Therefore we will consider the direct Ni-Ni interchain coupling irrelevant for a physical description of the material.

(ii) When Y, located between the chains, is completely substituted by a magnetic rare earth R, the system orders antiferromagnetically. The Néel temperature is typically smaller than the gap energy. For example, in the case R=Nd the $`T_N=48`$K and $`\mathrm{\Delta }11`$meV ($`127`$K). In the paramagnetic phase ($`T>T_\mathrm{N}`$), the Ni-chain excitations in R2BaNiO5 are virtually indistinguishable from those in the Yttrium compound at the same temperature, implying that very little should change in the model of the $`S=1`$ chains while the ordering should be driven by some kind of exchange Ni-R that induces an indirect coupling between the chains.

(iii) In the magnetically ordered phase, the staggered magnetization of the Ni sublattice has as $`T0`$ saturation values of $`1.01.6\mu _B`$ per ion in all the R2BaNiO5 studied to date . These values are clearly smaller than the classical result $`2\mu _B`$. This observation signals that the model describing these compounds has to properly retain the quantum fluctuations in the ordered phase.

(iv) Unlike Ni-R interactions, that are vital in order to induce a static magnetization on the Ni sites below $`T_N`$, direct R-R magnetic coupling can be disregarded in our model. Indeed, the ordering temperatures in R2BaNiO5 compounds are typically several tens of Kelvin. Direct coupling strength between rare earths in insulators is much weaker, usually of the order of one Kelvin. This is due to the fact that in rare earth species the magnetic $`f`$-electrons are strongly localized, what prevents efficient superexchange coupling. Unlike the case of Y2-xCaxBaNiO5 , where doping with Ca introduces hole carriers in the Ni-chains, this localized nature of f-electrons in the rare-earth keeps the S=1 chains free from charge carriers.

(v) In the crystal structure of R2BaNiO5, the site-symmetry for R3+ is very low. As a result, for each rare earth ion the degeneracy of its magnetic multiplet is lifted by crystal field effects. To investigate the low-energy part of the spectrum we only need to consider the lowest-energy orbital levels. In the case of the rare earth being a Kramers ion with half-integer total angular momentum, we shall model the magnetic rare earths by effective $`s=\frac{1}{2}`$ (pseudo)spins. In order to decide about the type of coupling between the Ni and the R sublattices, we rely on the experimental observation of the absence of dispersion in the crystal field excitations associated with R3+ ions. This fact suggests that the coupling between the R and Ni sublattices is very anisotropic and can be approximated by an Ising-type term in the Hamiltonian.

With all the ingredients above, we propose now the following Hamiltonian:

$`H=J{\displaystyle \underset{ij}{}}𝐒_{i,2j}𝐒_{i+1,2j}+J_c{\displaystyle \underset{ij}{}}S_{i,2j}^z(s_{i,2j1}^z+s_{i,2j+1}^z)`$ (1)
with J and Jc both positive. The first term in $`H`$ is a $`S=1`$ Heisenberg model along the Ni chains, the second is an Ising-like coupling with strength Jc between the $`S=1`$ Ni ions and the $`s=\frac{1}{2}`$ magnetic moment in the rare earth. The index $`i`$ runs along the chain direction and $`j`$ in the direction perpendicular to the chains. A sketch of the model is shown in Fig. (1). When $`J_c=0`$, $`H`$ should reproduce the physics of the Y2BaNiO5, i.e. basically independent $`S=1`$ chains.

We have studied this model numerically with the Quantum Monte Carlo (QMC) method using a quantum cluster algorithm, the Loop Algorithm . This method allows for an efficient sampling of the configuration space giving a specific prescription of how to perform global updates in quantum systems. Such updates involve clusters of spins with a size of the correlation length. The hamiltonian we studied is not frustrated and does not show the sign problem. All the Boltzmann weights can be taken positive after the conventional rotation around the z-axis of all the spins in one of the two sublattices of the S=1 system. The difficulty in reaching extremely low temperatures is a drawback of the QMC method when one is interested in characterizing exotic ground states, mostly near quantum critical points or strongly frustrated systems. None of these cases affects our calculations since our main goal is to compare with experimental results at temperatures, as we will see, perfectly accesible for QMC. We have measured various magnitudes that can be directly compared with the available experimental data and which, in some cases, are not accesible for other approaches like the MF/RPA approximation or the zero temperature Density Matrix Renormalization Group (DMRG) calculations. These magnitudes are: (i) Magnetizations in both sublattices (ii) Néel temperatures (iii) Correlation lengths and gaps in both subsystems.

Before discussing the results, we present a brief technical description of the numerical work and data treatment. QMC simulations are performed in finite lattices therefore an exhaustive finite size analysis of the results is necessary to compare with real experiments. We used periodic boundary conditions. The magnetization ($`M`$) was computed by directly measuring $`|M|`$ in order to avoid the averaging between configurations around the two ground states of the model in the ordered phase, and by using the relation $`M=lim_L\mathrm{}|M|`$. We considered for this extrapolation L$`\times `$L lattices of size L=8,16,24,32,48,64 spins. The Néel temperature was then determined by using the Binder parameter which is the fourth cumulant of the order parameter distribution. At $`T_N`$ this cumulant is independent of the size of the system, apart from subdominant corrections. We have observed that for values of $`J_c`$ smaller than $`0.1J`$ these corrections are due to the anisotropy in the scaling (i.e. different number of spins in $`a`$ and $`c`$ directions). We have checked this point by studying systems with different sizes $`L_a=32,48`$ spins and $`L_c=24,32,48`$ spins. For values of $`J_c0.3`$ and the precision needed for comparison with the experiment, $`T_N`$ can be computed by extrapolation from isotropic lattices. The statistical error bars are always smaller than the symbol sizes in all the figures presented. We measured all the magnitudes in $`10^5`$ Monte Carlo steps after thermalization. We have used units $`\mu _B=k_B=\mathrm{}=1`$.

The staggered magnetizations obtained from the QMC computations for both the Nickel ($`M_{\mathrm{Ni}}`$) and rare earth ($`M_\mathrm{R}`$) sublattices are in good agreement with the experimental results for Nd2BaNiO5  for a suitable choice of the model parameters in the Hamiltonian. For instance, for a value of the transverse coupling $`J_c=0.31J`$, we have that $`M_{\mathrm{Ni}}(T0)=0.79`$ that corresponds to the value 1.6$`\mu _B`$ observed experimentally. The Néel temperature for this value of $`J_c`$ is $`T_N=0.163J`$. We can take the gap in a single $`S=1`$ chain, $`\mathrm{\Delta }=0.410J`$ and compute the ratio between the two main energy scales in the system $`r=\frac{\mathrm{\Delta }}{T_N}`$. We observe that $`r`$ is very similar in the material ($`r`$=2.59) and in the model ($`r`$=2.51). In Fig. (2) we present the temperature dependence of the staggered magnetizations in the $`S=1`$ (open squares) and the $`s=\frac{1}{2}`$ (open circles) sublattices, $`M_{\mathrm{Ni}}`$ and $`M_\mathrm{R}`$ respectively, for the selected value of $`J_c`$ obtained from the QMC simulations. The magnetization results have been rescaled with respect to the maximum saturation values, i.e. $`1`$ for the $`S=1`$ system and $`0.5`$ for the $`s=\frac{1}{2}`$ system. The results can be compared with the experimental magnetic moment in the Ni (filled squares) and Nd systems (filled circles). Note the good agreement between the theoretical and experimental results at low temperatures. We also observe that the staggered magnetization in the $`S=1`$ sublattice is nearly temperature independent for $`T\frac{T_N}{2}`$. At T=40K, the Nd2NiBaO5 undergoes a magnetic structural transition  and therefore we expect only qualitative agreement with our model at temperatures near $`T_N`$.

Both sublattices behave qualitatively different in the ordered phase. While in the Ni sublattice the reduction of the staggered magnetization by purely quantum fluctuations is important, the R sublattice is fully saturated to the classical value $`M_\mathrm{R}=0.5`$ at T=0.

In order to understand the origin of the ordered phase of this quantum model we have studied the relation between the magnetizations in both sublattices for different values of $`J_c`$. In Fig. (3)(a) we show a plot of the QMC results at $`T<T_N`$ for $`M_R`$ as a function of $`M_{Ni}`$ for various values of $`J_c`$. We observe that the behavior of the magnetizations is well reproduced by the function:

$`M_\mathrm{R}=M_0\mathrm{tanh}(\alpha \beta M_{\mathrm{Ni}})`$ (2)
where $`\beta =1/k_BT`$, $`M_0`$ is the effective moment of the rare earth ion and $`\alpha `$ is a linear function of $`J_c`$ in the range $`0.09<J_c<0.5J`$ as observed from our QMC results. Eq. (2) can be obtained by a mean-field approach as shown by Zheludev et al.. These authors considered that the behavior of the R2BaNiO5 compounds in the ordered phase could be described in terms of a $`S=1`$ chain of Ni ions in a staggered magnetic field induced by the $`s=\frac{1}{2}`$ magnetic rare earth ions and, reciprocally, the otherwise free spins $`s=\frac{1}{2}`$ see the mean field produced by the neighboring $`S=1`$ chains.

The excellent agreement between this mean-field approach and our QMC results indicates that our microscopic model also fulfills the property of viewing the effect of the staggered magnetization on the rare earth ions as an effective magnetic field on the Ni subsystem and viceversa.

In Fig. (3)(b) we show the relation between $`M_{\mathrm{Ni}}`$ and the staggered magnetic field induced by $`M_\mathrm{R}`$ obtained by our QMC calculations for various T values. The extrapolation $`T0`$ reproduces the DMRG results for one $`S=1`$ chain in a staggered magnetic field at $`T=0`$ obtained by Yu et al.. Since the R sublattice is fully polarized at $`T=0`$ no rescaling is necessary to compare our data with the DMRG results. The quality of the fits in Fig. (3)(a) and the data collapse in Fig. (3)(b) suggests a technical procedure to extrapolate QMC data down to zero temperature. This procedure is obvious for magnetizations but could be extended for other thermodynamic magnitudes for Hamiltonians similar to the one proposed here.

In the context of the effective field interpretation, an analysis of the temperature dependence of $`M_\mathrm{R}`$ obtained with our QMC simulations indicates that though $`M_\mathrm{R}`$ saturates to the classical value at T=0 (see Fig. (2)), as expected if we consider this subsystem as a sublattice of free spins in the presence of an external field, this behavior cannot be reproduced by a purely classical mean field model in the R sublattice, namely mean field just proportional to the magnetization in the R sublattice and independent of the Ni magnetization. This fact is signaling that quantum fluctuations in the Ni sublattice are essential to describe the magnetization in the R sublattice.

We discuss now the spatial spin-spin correlations in the paramagnetic phase. In Fig. (4) we plot the inverse of the correlation length $`\xi ^1`$ along the Ni chains (circles) and along the R chains parallel to the Ni chains (squares) for $`T<\frac{\mathrm{\Delta }}{2}`$. For comparison, we present $`\xi ^1`$ in a single $`S=1`$ spin chain (triangles). In this range of temperatures it is well known that the correlation length in the $`S=1`$ chain is $`\xi =\frac{\mathrm{}c}{\mathrm{\Delta }}`$ where $`c`$ is the spin-excitation velocity. The correlation length is computed in all the cases by fitting the z-component of the spin-spin correlation function to the law:

$`|S_0^zS_l^z|=A\mathrm{exp}({\displaystyle \frac{l}{\xi }})l^\eta `$ (3)
($`A`$ and $`\eta `$ are fitting parameters) directly in a lattice of 64$`\times `$64 spins and for $`J_c=0.18`$ to keep $`T_N`$ and $`\mathrm{\Delta }`$ well separated in energy. Near $`T_N`$, the critical modes become gapless and we observe how the gap closes in both subsystems following a law $`\xi ^1=K(TT_N)^{\frac{1}{2}}`$. K is the only free parameter in the fit presented in Fig. (4) as a solid line, since $`T_N`$ is computed independently. As the temperature increases, the correlation length in the Ni subsystem approaches the correlation length of a single chain.

The behavior of $`\eta `$ is completely different in both sublattices. While in the R sublattice the correlation function $`S_0^zS_l^z`$ can be fitted with $`\eta =0`$ for all the temperatures showed in Fig. (4), in the Ni sublattice $`\eta `$ approaches the value $`0.5`$ in Eq. (3) as the temperature is reduced in the paramagnetic phase. We have analized this behavior and observed that the temperature dependence of $`\eta (T)`$ in the spin-spin correlation function Eq. (3) is very similar to that of $`\eta `$ in a single chain in which Eq. (3) shows a slow cross-over from $`\eta =\frac{1}{2}`$ at T=0, more precisely $`S_0^zS_l^z`$ is then described by the modified Bessel function $`K_0`$ , to $`\eta =0`$ when $`T>\mathrm{\Delta }`$ . The temperature dependence of $`\eta (T)`$ is shown in Fig. (5).

The picture arising from the study of the spin-spin correlation functions is that the modes becoming gapless at $`T_N`$ can be described at a mean field level. Moreover, since both the value of the correlation length and the power law corrections to the exponential behavior are very similar for the Ni subsytem and for an independent $`S=1`$ chain, we can conclude that the modes remaining gapped in the system are very much like the conventional Haldane excitations.

In summary, we have shown that QMC calculations based on the proposed microscopic model are able to reproduce the observed behavior of R2BaNiO5 mixed-spin quantum antiferromagnets remarkably well. They provide a solid numerical basis for the MF/RPA model, and make additional predictions regarding the magnetic nature of both spin species that go beyond the simplistic picture of a MF/RPA model. We hope that the proposed approach will be particularly useful in the study of dilute (RxY1-x)2BaNiO5 compounds in the limit of small $`x`$, where the MF model becomes totally inadequate.

It is a pleasure for us to acknowledge discussions with C. Gros, F. Mila and H. Rieger. Oak Ridge National Laboratory is managed by UT-Battelle, LLC for the U.S. Department of Energy under contract DE-AC05-00OR22725.

\*Previous address: Physics Department , Brookhaven National Laboratory, Upton, NY 11973-5000, USA.


text:
# 1 Introduction

## 1 Introduction

Nowadays, it has been a matter of increasing interest  to apply the fundamentals of the theories of complex systems in many different disciplines, not only in physical sciences, but even in social sciences, from economy to education  or sociology. The main point is that social systems, like natural ones, are constituted of great number of individuals, which - generally - have local interactions between them. Sometimes, social networks behaviour can be determined also by the action of external actors, which might be mimicked by external fields in our model.

Elections are processes where many individuals interact between them. It is a dynamical convincing process, where we have at the same time the interaction between neighbours and external influence (political advertising, campaigns etc). In Brazil, in proportional elections (deputies or city councillors) the voters vote directly for the candidates and not for the parties. They can vote for a party, but it is not frequent. Some elections occur with a large number of voters: In some states or in the largest cities one has a number of voters in the order of magnitudes of millions or tens of millions. So, these elections are a social phenomenon which presents the basic characteristics of complex systems. One of these features is that they are scale-free phenomena. This feature has been observed by Costa-Filho et al. , who showed that the distribution of the number of votes obtained by different candidates for the 1998 elections in Brazil follow a power law distribution, with exponent $`1.0`$. The same result can be obtained for the whole country or for other proportional elections, which shows that this is a very robust result.

Taking into account that elections are processes where a vote is supposed to be obtained as a result of convincing arguments, it can be compared with a physical process of clustering. However, contrary to the results obtained by Costa-Filho et al, usual models of formation of clusters (as percolation, for instance ) may give exponents $`2.0`$ (square lattice) for the numbers of clusters as a function of the cluster size.

Recently, we have introduced  a model for proportional elections. Our model is based on the Sznajd model proposed to simulate the process of formation of opinion. However, different from other models , where the influence flows inward from the border to the center (like in the majority games, where the site in the middle takes the state of the majority of neighbouring sites), in the Sznajd model  one has an outward flow of influence. It thus differs from e.g. bootstrap percolation  or other cellular automata (for a computational review see ) where the site in the center behaves according to a rule determined by its neighbours. Nevertheless, the dynamics of the Sznajd model is quite similar (except for isolated sites) to that of spinodal decomposition of the Ising model  at low temperatures, as was shown in : Starting from a random distribution, large domains form where nearly all sites have the same state. Finally, one domain will cover the whole lattice. Thus to get the desired results we will look at intermediate times when there are still many different domains or correlated sites in the system.

In the Sznajd model, small sets of people influence the opinions of their nearest neighbours if and only if all people within the original set agree. On a chain, this set is a bond with two people at its ends . On the square lattice with  or without  disorder, it can be such a bond or a plaquette of four neighbouring spins (people). This plaquette rule is called rule Ia in  and the bond rule, used in the present paper, is called rule IIa in . Thus if all four plaquette members or two bond members share the same opinion, they convince their neighbours of this opinion; otherwise the neighbours remain unchanged.

However, as shown by Barabasi and Albert recently , social relations must be represented by networks instead of lattices. Networks of interactions (www, author’s collaboration in scientific papers, actor’s collaborations in films) show the common feature of scale-free behaviour. In order to represent this main feature, Barabasi and Albert introduced a model for evolving networks. Numerous papers used this model for a variety of purposes, e.g. . Starting with few nodes (which may represent actors, authors, web sites) connected to each other, more and more nodes are added to the network, each node connecting to an already connected one, with the probability to connect to a node being proportional to the number of previous nodes which are already connected to it.

In Brazil the voting process is much more based on the relation between candidate/voter than on the parties. Thus, our idealized version of the voting process can omit the role of the parties. Another aspect is that it was clear for us that the candidates do not start with the same social weight. This determines the result of the elections, since candidates with more social visibility or better conditions to campaign are more likely to be elected. So, we have to introduce some differentiation between the candidates in the beginning of the simulation.

In this paper, we performed simulations on a three-dimensional version of the previous one-dimensional  and two-dimensional  Sznajd model. We combine it with a network model for elections based on the model of Bernardes et al . Unlike this first version  and its three-dimensional variant where a probability to convince had to be introduced, eq.(1) below, in order to produce some differentiation between the candidates, with the Barabasi network the same result is obtained from the combination of the different number of neighbours of the nodes without this probability.

In the next section we present the models we have used, followed by the results. Both of these models use networks connecting the voters; one network is a simple cubic lattice, the other a Barabasi network. After that, we conclude.

## 2 Models and Results

We have simulated two models in the present work. The first one is a 3d version of that simulated previously . The second is a Barabasi network version.

### 2.1 Simple cubic lattice

In this work, we used a modified version of the Sznajd model (rule IIa in ): A pair of neighbours in agreement convinces its ten nearest neighbours to the same opinion. A cubic lattice of size $`L\times L\times L`$ represents the set of voters. A number $`N_{tot}`$ of candidates, $`N_{tot}L^3`$, is fixed in the beginning of the simulation. The value $`n=1,2,\mathrm{},N_{tot}`$ of a site $`S`$ on the lattice will represent that this voter prefers that candidate $`n`$. The model has two different stages: First, we produce the initial condition and, after that, we perform the simulation of the electoral campaign (only voters can influence other voters, a la Sznajd). As in real elections, we do not wait for a kind of equilibrium state, but count the votes at some intermediate time. Basically what we are doing is the analysis during the transient time. As in real elections, the candidates have different initial chances of being voted for (representing more money for electoral campaigns, more initial visibility etc.). This is modelled by a probability $`P_c`$ of convincing, calculated from the label $`n`$ of the candidate

$$P_c=(n/N_{tot})^2.$$
$`(1)`$
It means that the higher is the label $`n`$ of a candidate, the higher is the probability of convincing a voter.

In the first stage, we started with all the sites with value zero, meaning that there are no committed voters. Then, we visit all the sites exactly once, in random order. For each visit, we try to convince the voter to adopt a candidate, chosen at random. A random number $`r`$ is generated and compared with $`P_c`$. If $`rP_c`$ the candidate is accepted by that voter. If the candidate convinces the voter, this voter tries to convince the neighbouring sites. Once again, we throw the dice and compare a new random number with $`P_c`$. If successful, $`rP_c`$, the voter will try to convince the neighbourhood as follows: We check all the six neighbouring sites; for each that has the same value of the candidate chosen before, all the ten neighbouring sites of this bond of two sites will assume the same value (as in the usual Sznajd prescription). If nobody has chosen the same candidate, only the originally selected voter is committed to this candidate.

In the second stage, a usual Sznajd process is performed without using the complication from the probability $`P_c`$. (We thus assume all voters to be equal and restrict the probability $`P_c`$ to describe the convincing power of the candidates only.) We go to random sites on the lattice. A neighbouring site is chosen at random and we check if the two sites have the same value (they prefer the same candidate). In that case, all the ten neighbours change to vote in that candidate.

Figure 1 shows, just as in real life , deviations from a simple power law for both very large numbers of votes and very small numbers. In between, however, the simulations are compatible with the hyperbolic law

$$N(v)1/v$$
$`(2)`$
observed in reality for the number $`N`$ of candidates having $`v`$ votes each. (Here and in Fig.2 below the bin size for $`v`$ increases by a factor 2 for each consecutive bin.)

\[In , for the square lattice and assumption (1), two exponents are fitted onto the data: One for the distribution after the first stage, and one for the distribution during the second change. When assumption (1) is generalized from $`n^2`$ to $`n^x`$, then the first exponent depends strongly on $`x`$ (which can be explained by a simple analytical relation between this exponent and $`x`$), while the second exponent depends much less on $`x`$.\]

### 2.2 Network version

In this model, we first create a network of interacting nodes by using the basic Barabasi-Albert prescription. We fix an initial number of nodes, each one connected to the others. In the present work, the minimum number of connections of a node is $`m=5`$. So, in the beginning of the simulation we have 6 nodes, in order that each one can be connected to the 5 others. After that, more and more nodes are added to the network. A new node has a probability to be connected to a previous node proportional to the number of nodes that are already connected to this previous node. Thus the growth probability at any existing node is proportional to the number of nodes already connected to it. We will no longer need assumption (1) and have replaced this assumption and the regular lattice by the Barabasi network without such an assumption.

After preparing the network, we start with the election process, which is now different from that in the previous three-dimensional model. The first step is the distribution of candidates. Again, the state of a node, that means, the value $`n`$ of a node on the network, represents that this voter has given the preference to that candidate $`n`$. Thousand candidates are distributed at random, disregarding the number of connections of a node, i.e. we pick a node at random from the half million nodes to which we let the network grow, and then a candidate at random. Now, the campaign starts. At each time step we visit all the nodes. For each node, we have the following process:

* If a node $`i`$ has already selected preference for a candidate, we choose a connected node $`j`$ at random. If node $`i`$ has no candidate ($`n=0`$), we go to another randomly selected node.
* If node $`j`$ has the same candidate as node $`i`$, they try to convince all the nodes connected with them. The probability to convince others for each of the two nodes is now inversely proportional to the time-independent numbers of nodes connected with it, meaning that each node convinces - on average - one other node at each process.
* If node $`j`$ has no candidate, node $`i`$ tries to convince it to accept its own candidate, with the same probability as described above.
* If node $`j`$ has a different candidate from node $`i`$, we skip to another node $`i`$.

Again, as described above for the 3d version, we do not wait for a equilibrium state. It is important to mention that, different from a square lattice, where an equilibrium state is reached in a time proportional to the number of sites, in the network it is reached rapidly, after about $`10^2`$ iterations. In both cases, in final equilibrium all the sites have the same state.

Figure 2 shows that again except for the smallest and the largest numbers $`v`$ of votes, the hyperbolic law (2) is obeyed well at intermediate times $`t=40`$.

## 3 Summary

Whether we simulate the election process on a square lattice, a simple cubic lattice or a Barabasi network, we recover the same hyperbolic law as found in real elections. Our simulations on a Barabasi network have the advantage that we no longer need assumption (1) for the purpose of getting a realistic vote distribution with decay exponent 1 in the center.

Our study has shown that the hyperbolic law observed empirically is a rather robust consequence of our modifid Sznajd model, since we found it first on the square lattice  and now on both the cubic lattice and the Barabasi network. Either we use a regular lattice and assumption (1), or we use the Barabasi network without assumption (1); the final results are similar. The fact that the hyperbolic law is observable on the Barabasi network, which is a more realistic model of social interactions than the lattices, provides evidence that the Sznajd model may well capture important aspects of the voting mechanism. The advantage of using the Barabasi network instead of regular lattices is not only that it is more realistic but also that we can drop assumption (1) which is a kind of fine tuning the system to criticality. Of course, the Barabasi model, and the related assumptions are also rules (as we have rules when constructing a lattice too) but the difference is similar to what we have for usual and self-organized criticality: Assumption (1) puts in some exponent at the beginning, through $`(n/N_{tot})^2`$, on which the final exponent depends somewhat (see end of section 2.1), while the Barabasi growth process leads by itself to a power-law distribution of the number of connections for a node, and combined with the Sznajd model gives the desired vote distribution with its intermediate power law.

Moreover, rule (1) was introduced ad hoc to explain the election results, while the rules of  were stated before, independent of the present application. We are not aware of other voter models explaining the hyperbolic law found empirically in .

Acknowledgments: We thank Ana Proykova for suggesting the 3d simulation. ATB acknowledges the hospitality of the Institute for Theoretical Physics from the University of Cologne. This work was partially supported by the Brazilian Agencies CNPq, FINEP and by the Hungarian OTKA T029985.