text:
# The Interstellar Rubidium Isotope Ratio toward Rho Ophiuchi A

## 1 Introduction

Studies on rubidium in stellar atmospheres and the interstellar medium (ISM) are few, but it is not for a lack of interest in its abundance. The production of Rb involves the neutron capture $`s`$\- and $`r`$-processes. Rubidium production by the $`s`$-process occurs through the ‘weak’ process in the He- and C-burning layers of massive stars and through the ‘main’ process in the He-shell of low-mass AGB stars. Supernovae from core collapse of massive stars are the likely source of the $`r`$-process. Analysis of the Solar System abundances provides estimates of the fractional responsibility of these neutron capture processes for the Rb isotopes 85Rb and 87Rb. Our estimates drawn from averaging the results of Beer, Walter, & Käppeler (1992), Raiteri et al. (1993), and Arlandini et al. (1999) are that 85Rb is about 35% $`s`$\- and 65% $`r`$-process in origin with the weak $`s`$-process being one third of the $`s`$-process contribution. For 87Rb, the fractions are about 70% from the $`s`$\- and 30% from the $`r`$-process with a minor contribution from the weak $`s`$-process. As a result of the different relative contributions of the $`s`$\- and $`r`$-processes to the two isotopes, a measurement of the Rb isotope ratio is a clue to the history of heavy element nucleosynthesis.

Unfortunately, rubidium is a difficult element to measure. It is potentially measureable only in cool stars. Resonance lines of Rb I are detectable at 7800 and 7947 Å, but blending with atomic and/or molecular lines adds an unwelcome difficulty to an abundance analysis. Rubidium abundances for unevolved stars were reported by Gratton & Sneden (1994) and Tomkin & Lambert (1989). Measurements of the abundance for giant stars were obtained by Lambert et al. (1995), Abia & Wallerstein (1998), and Abia et al. (2001). The primary aim of these latter studies was to exploit the use of Rb as a neutron densitometer for the $`s`$-process site (Lambert et al. 1995). This use arises from the role of 85Kr as a branch in the $`s`$-process path. The branch directly affects the isotope ratio, 85Rb/87Rb, but it is not directly measureable from even very high resolution stellar spectra (Lambert & Luck 1976) because the lines are too broad. Earlier interstellar measurements yielded upper limits (Federman et al. 1985) because the Rb abundance is quite small ($`2.5\times 10^{10}`$ in meteorites; Anders & Grevesse 1989), and it is derived from Rb I, not the dominant form Rb II. The only interstellar detection so far comes from observations on two heavily reddened stars in Cygnus (Gredel, Black, & Yan 2001). \[An earlier reported detection by Jura & Smith (1981) could not be confirmed by Federman et al. (1985).\] Here, we present high-resolution spectra revealing not only another interstellar detection, but also the first measurement determining the 85Rb/87Rb isotope ratio for extrasolar gas.

## 2 Observations and Analysis

The star $`\rho `$ Oph A is an ideal target for a study on Rb isotopes. It is relatively bright ($`V`$ $`=`$ 5.0), moderately reddened with $`E`$($`BV`$) $`=`$ 0.47, and has one main interstellar component (e.g., Lemoine et al. 1993). The weak lines of Li I (Lemoine et al. 1993) and K I $`\lambda `$4044 (Crutcher 1978) toward this star also show more absorption than is typical for a diffuse cloud.

The data on Rb I $`\lambda `$7800 were acquired with the 2dcoudé spectrograph (Tull et al. 1995) on the Harlan J. Smith 2.7 m telescope at McDonald Observatory in 2003 May. We used the high-resolution mode with echelle grating E2, centered on 7444 Å in order 46. The spectra were imaged onto a $`2048\times 2048`$ Tektronics CCD (TK3). The 145 $`\mu `$m slit provided a resolving power of 175,000 as determined from widths of lines in the Th-Ar comparison spectrum, sufficient to discern broadening of the interstellar lines. We obtained calibration exposures for dark current the first night of the run, bias correction and flat fielding each night, and Th-Ar spectra every 2 to 3 hours during the night. Our primary target, $`\rho `$ Oph A, was observed for a total of 8h, with individual exposures of 30m per frame; this procedure minimized the deleterious effects caused by cosmic rays. In addition, the unreddened star $`\alpha `$ Vir was observed for a total of 100m centered on the slit to check for contamination from weak telluric lines and CCD blemishes not removed during the flat-fielding process. No contamination is present. An additional 60m was utilized to trail $`\alpha `$ Vir along the slit to act as a stellar flat field. The analysis described below utilized spectra of the K I line at 4044 Å that were acquired for another project (Knauth, Federman, & Lambert 2004). While details will be presented in Knauth et al. (2004), we note that the instrumentation and observing procedure were basically the same as those for Rb I, except echelle grating E1 was used. There is at most a 0.1 km s-1 difference in line velocities between the two setups according to the dispersion solutions found from Th-Ar spectra.

Standard routines within the IRAF environment were used to extract one-dimensional spectra that were Doppler-corrected and normalized to unity. The Rb I spectrum for the interstellar gas toward $`\rho `$ Oph A is displayed in Fig. 1. The appearance of two ‘features’ arises from the 62 mÅ hyperfine splitting in 85Rb combined with the stronger (redder) hyperfine component in 87Rb. The rms deviations in the stellar continuum yield a signal-to-noise ratio per pixel of 1200 for $`\rho `$ Oph A; there are 2.9 pixels per resolution element.

As in our earlier work on the Li isotope ratio acquired at the same resolution (Knauth et al. 2000), the Rb I and K I lines were fitted to extract column densities for neutral 85Rb, 87Rb, and K. The relevant atomic data used for input (Morton 1991, 2000) and the resulting equivalent widths ($`W_\lambda `$) and column densities ($`N`$) are given in Table 1. There is one main component at $`V_{LSR}`$ of 1.9 km s-1 (Lemoine et al. 1993; Pan et al. 2003), which at ultra-high-resolution is seen as two in K I $`\lambda `$7699 spectra (Welty & Hobbs 2001). Another component at 3.5 km s-1 appears in K I $`\lambda `$7699 (Welty & Hobbs 2001; Pan et al. 2003), having $``$ 20% of the column of the main component. This second component is marginally seen (at the 2-$`\sigma `$ level) in our K I $`\lambda `$4044 spectrum, but there is no evidence for it in the Rb I line. The uncertainties in column densities were inferred from a map of chi-squared confidence levels (Fig. 2). While optical depth effects are not important for these very weak lines, the profile fitting code (see Knauth et al. 2003) determines the $`b`$-value as well (from the difference between measured line width and the instrumental width determined from Th-Ar lines). We obtained $`b`$-values of 1.0 and 0.8 km s-1 for Rb I and K I; the slight difference is not significant at our spectral resolution.

## 3 Results

The primary result of our study is the determination of a Rb isotope ratio for the main cloud toward $`\rho `$ Oph A. The best fit (with a reduced chi squared of 1.43 for 47 degrees of freedom) yields a 85Rb/87Rb ratio of $`1.21\pm 0.30`$ – see top panel of Fig. 1. This differs significantly from the meteoritic value of 2.59 (Anders & Grevesse 1989). We attempted to fit our data with the Solar System ratio as well (bottom panel of Fig. 1). The reduced chi squared is worse (2.08) and an F-test shows that there is 50% confidence in this fit. Examination of the lower panel reveals where the difference lies: Neither the blue nor red shoulders of the lines are fit well using a Solar System ratio. Furthermore, the goodness of fit can be judged by variations in the residuals (data minus fit) outside and inside the absorption profile. The variation in the lower panel is not as consistent across the spectrum. Other syntheses with the Solar System ratio were attempted without success. Because both shoulders show more absorption than is predicted by this ratio, addition of the redder 3.5 km s-1 component would not improve the fit. Broader lines could fill in the shoulders of the observed Rb I profile, as could a two-component fit upon shifting the profile by about 1 resolution element. \[A $`b`$-value as large as 1.8 km s-1 is possible from our $`\chi ^2`$ analysis, though it is not consistent with other data (e.g., Welty & Hobbs 2001).\] However, the 85Rb hyperfine components are broadened as well, yielding a chi squared that is almost double the $`\chi ^2`$ of the best fit in both cases. Finally, we note that the lines representing different values of 85Rb/87Rb in Fig. 2 favor a ratio between 1.0 and 1.5. The Solar System ratio of 2.6 differs by 2- to 3-$`\sigma `$ in each quantity (column of Rb isotope), consistent with the 4.5-$`\sigma `$ difference seen in the uncertainty listed for the 85Rb/87Rb ratio.

The Rb/K elemental ratio provides a means to determine for this interstellar cloud whether the 85Rb abundance is lower or the 87Rb abundance is higher relative to the abundances in meteorites. Since Rb I and K I are minor interstellar species, ionization balance is required to extract an elemental abundance. The elemental ratio does not depend on electron density, and we assume that the depletion levels for alkalis are the same (see Welty & Hobbs 2001; Knauth et al. 2003). Then for the Rb/K ratio, we have

$$\frac{A_g(\mathrm{Rb})}{A_g(\mathrm{K})}=\left[\frac{N(\mathrm{Rb}\mathrm{I})}{N(\mathrm{K}\mathrm{I})}\right]\left[\frac{G(\mathrm{Rb}\mathrm{I})}{G(\mathrm{K}\mathrm{I})}\right]\left[\frac{\alpha (\mathrm{K}\mathrm{I})}{\alpha (\mathrm{Rb}\mathrm{I})}\right],$$
(1)
where $`G`$(X) is the photoionization rate corrected for grain attentuation and $`\alpha `$(X) is the rate coefficient for radiative recombination. Comparison of the theoretical calculations for radiative recombination of Rb and K (Wane 1985; Péquignot & Aldrovandi 1986; Wane & Aymar 1987) reveals that the coefficients are the same to within about 10%, confirming the assumption made in other interstellar studies (Jura & Smith 1981; Federman et al. 1985; Gredel et al. 2001). The photoionization rates for Rb I and K I are $`3.42\times 10^{12}`$ and $`8.67\times 10^{12}`$ s-1, respectively. The cross sections for Rb I are from the theoretical calculations of Weisheit (1972) for 1150 to 1250 Å, scaled to the measurements of Marr & Creek (1968) at longer wavelengths, while those for K I are from the measurements of Hudson & Carter (1965, 1967), Marr & Creek (1968), and Sandner et al. (1981). Since the ionization potentials for Rb I and K I are similar, 4.18 and 4.34 eV, there is no appreciable differential attenuation. We use this fact to account for Rb I absorption below 1150 Å; based on the total rate for K I, we applied a 13% ‘correction’ to $`G`$(Rb I). The result is a Rb/K ratio of $`(1.3\pm 0.3)\times 10^3`$ compared to the Solar System ratio of $`(1.9\pm 0.2)\times 10^3`$ (Anders & Grevesse 1989).

Limits or values of the Rb/K ratio are available from earlier studies of interstellar Rb as well. For the clouds toward $`o`$ Per, $`\zeta `$ Per, and $`\zeta `$ Oph, the most conservative limits on Rb I absorption (Federman et al. 1985) are used with our more recent determinations of the K I column density from the weak line at 4044 Å (Knauth et al. 2000, 2003). The 3-$`\sigma `$ upper limits for $`N`$(Rb I) are $`2.8\times 10^9`$, $`1.9\times 10^9`$, and $`3.7\times 10^9`$ cm-2, respectively. The limits on the Rb/K ratio then become $`1.2\times 10^3`$, $`1.0\times 10^3`$, and $`1.8\times 10^3`$. The Rb I detections toward Cyg OB2 Nos. 5 and 12 (Gredel et al. 2001) can be combined with results from high-resolution K I $`\lambda `$7699 spectra (McCall et al. 2002). For the three main molecular components ($`+4.0`$, $`+6.5`$, and $`+12.3`$ km s-1), we derive values for $`N`$(K I) of $`9.4(7.7)\times 10^{11}`$, $`8.1(12)\times 10^{11}`$, and $`6.9(16)\times 10^{11}`$ cm-2 for the gas toward No. 5 (No. 12) when adopting a $`b`$-value of 1 km s-1. The respective Rb/K ratios are $`1.4\times 10^3`$ and $`1.2\times 10^3`$. (The uncertainties for these ratios are hard to quantify because the Rb I component structure is not known and because the K I column densities from $`\lambda `$7699 are very susceptible to small changes in adopted $`b`$-value.) These comparisons suggest that the interstellar Rb/K ratio may be lower than the meteoritic abundance throughout the solar neighborhood.

## 4 Discussion

Insight into the nucleosynthesis of Rb is obtainable from the ratios 85Rb/K and 87Rb/K for the main 1.9 km s-1 component toward $`\rho `$ Oph A. We find 85Rb/K = $`(0.73\pm 0.12)\times 10^3`$ and 87Rb/K = $`(0.60\pm 0.13)\times 10^3`$, assuming that $`G`$(Rb I) and $`\alpha `$(Rb I) are not dependent on isotope. The Solar System ratios are 85Rb/K = $`(1.37\pm 0.14)\times 10^3`$ and 87Rb/K = $`(0.53\pm 0.06)\times 10^3`$. These ratios show that 85Rb is underabundant in the gas by a factor of about two, but 87Rb has about the Solar System abundance. In other words, the 85Rb/87Rb and Rb/K ratios toward $`\rho `$ Oph A are low because there is less 85Rb. The ‘missing’ 85Rb is comparable to predictions for the $`r`$-process component in stellar models (Beer et al. 1992; Raiteri et al. 1993), $``$ 65%. Moreover, the models (Beer et al. 1992; Raiteri et al. 1993; Arlandini et al. 1999) indicate that 87Rb arises mainly from the $`s`$-process. We are led to believe that the lower interstellar abundance for 85Rb is due to less $`r`$-process synthesis.

It is widely assumed that the site of the $`r`$-process is the very deep interior of a Type II supernova. Potassium is also a product of a SN II; it is synthesized in explosive oxygen burning (e.g., Clayton 2003). These different sites and quite different synthesis mechanisms make it likely that the $`r`$-process yield of a Rb isotope relative to the yield of K can vary; it may depend for example on the mass, metallicity, and rotation of a SN II’s progenitor.

The proposal that local interstellar gas may be deficient in $`r`$-process products relative to Solar System material is curiously supported by interstellar abundances of other lightly depleted or undepleted elements with a $`r`$-process contribution. We emphasize directions probing material toward Sco OB2, which includes $`\rho `$ Oph A. For Solar System material, about 50% of Kr, 50% of Cd, and 30% of Sn are attributable to the $`r`$-process (Beer et al. 1992; Raiteri et al. 1993; Arlandini et al. 1999). The interstellar Kr abundance within a few hundred parsecs of the Sun is constant at 60% of the solar-wind value (Cartledge et al. 2001, 2003). Though there are fewer determinations of the interstellar Cd abundance, it is constant at 80% of the Solar System value (Sofia, Meyer, & Cardelli 1999). Sofia et al. find a Sn abundance that is essentially solar, if not slightly greater than solar, for sight lines showing little depletion in other elements. Arsenic, which is mostly $`r`$-process (Beer et al. 1992; Raiteri et al. 1993) and has only 1 isotope, is seen currently only in three reddened directions (Cardelli et al. 1993; Federman et al. 2003). It has a higher condensation temperature (see Cardelli et al. 1993) and thus it has more depletion, especially in reddened directions. Measurements of As II in lightly reddened sight lines, where little depletion is expected, appear necessary to confirm our hypothesis.

Variations in space and time of the mix of $`s`$\- and $`r`$-process products in the Galaxy’s interstellar medium should be reflected in the abundances of the heavy elements in unevolved stars. Relative abundances of elements such as Sr, Y, and Zr and heavier elements such as Ba, Ce, and Eu for dwarfs of the thin disk do not vary from star to star by more than the measurement errors (Edvardsson et al. 1993; Chen et al. 2000; Reddy et al. 2003). These measurements would seem to exclude the larger variations expected for stars formed from gas with a large amplitude in the ratio of $`s`$\- to $`r`$-process products. Our interstellar results suggest that differences in $`n`$-capture nucleosynthesis may be local to the solar neighborhood.

## 5 Concluding Remarks

Our estimate of the Rb isotope and the Rb/K ratios combined with the interstellar abundances of Kr, Cd, and Sn raise the intriguing possibility that the $`r`$-process products for such light $`n`$-capture elements are relatively underrepresented in the local interstellar medium relative to their contribution in the 4.5 Gyr older Solar System material. While models of the $`s`$-process are becoming quite sophisticated and generally give similar results, differences limit the ability to make definitive statements. For these rarer elements, the Solar System abundance is sometimes a matter of debate; Kr is an example. The precision of the atomic data is another factor. For Sn II, the $`f`$-values adopted by Sofia et al. (1999) are now known to be valid (Schectman et al. 2000; Alonso-Medina, Colón, & Matińez 2003). A key test of our hypothesis is the measurement of accurate interstellar abundances for additional elements, particularly those derived mainly from the $`r`$-process.

We thank Dr. B. McCall, who kindly provided us with his K I spectra, and Dr. D. York for helpful suggestions. This work was supported in part by NASA LTSA grant NAG5-4957 to the Univ. of Toledo.


format: text

text:
# A First Principles Density-Functional Calculation of the Electronic and Vibrational Structure of the Key Melanin Monomers1footnote 11footnote 1This work relates to Department of the Navy Grants N00014-03-1-4115, N00014-02-1-1046 and N00014-03-1-4116 issued by the Office of Naval Research International field office. The United States Government has a royalty-free license throughout the world in all copyrightable material contained herein.

## I Introduction

The melanins are an important class of pigmentary macromolecule found throughout the biosphere Prota (1992). Pheomelanin (a cysteinyl-dopa derivative) and eumelanin (formed from 5-6-dihydroxyindolequinone and other indolequinones) are the predominant forms in humans, and act as the primary photoprotectant in our skin and eyes. Consistent with this role, all melanins show broad band monotonic absorption in the UV and visible in the range 1.5 to 5 eV Wolbarsht et al. (1981). In contrast, other biomolecules such as proteins and nucleic acids show only well defined absorption peaks around 280 nm (4.5 eV) and little absorption below that van Holde et al. (1998). They are also efficient free radical scavengers and antioxidants Prota (1992). In direct contradiction with these photoprotective properties, both pheomelanin and eumelanin are implicated in the development of melanoma skin cancer Hill (1995). For this reason, the photophysics, photochemistry and photobiology of melanins are subjects of intense scientific interest.

Despite work over several decades (with respect to eumelanin in particular), the more general structure-property-function relationships that control the behaviour of these important bio-macromolecules are still poorly understood Stark et al. (2003). The melanins are difficult molecules to study: they are chemically and photochemically stable, and are virtually insoluble in most common solvents. It is fairly well accepted that eumelanins are macromolecules of the various redox forms of 5,6-dihydroxyindolequinone (DHI or HQ) and 5,6-dihydroxyindole 2-carboxylic acid (DHICA) Prota (1992); Ito (1986). However, major questions still remain concerning their basic structural unit Zajac et al. (1994). Two opposing schools of thought exist: i) that eumelanins are composed of highly cross-linked extended hetero-polymers based upon the Raper-Mason scheme Prota (1992), and ii) that eumelanins are actually composed of much smaller oligomers condensed into 4 or 5 oligomer nano-aggregates Clancy et al. (2000). Clearly, this is a fundamental issue, and is the starting point for the construction of consistent structure-property-function relationships.

The answer to this question also has profound implications for our understanding of the condensed phase properties of melanins. In 1960 Longuet-Higgins Longuet-Higgins (1960), in a landmark paper, proposed that many of the physical properties of melanins could be understood if they were semiconductors. This proposition was lent further support by Pullman and Pullman Pullman and Pullman (1961) in 1964 who applied molecular orbital theory in its simplest form (the Hückel approximation) to 5,6-dihydroxyindole and one particular dimer type. This theoretical work was followed in 1974 by the experimental observations of McGinness, Corry and Proctor who demonstrated that a pellet of melanin could be made to behave as an amorphous electrical switch McGinness et al. (1974). They postulated that these materials may indeed be disordered organic semiconductors consistent with the Longuet-Higgins theory, and the newly developed models of amorphous inorganic semiconductors Mott and Davis (1971). Several studies since have also claimed to show that melanins in the condensed solid-state are semiconductors Crippa et al. (1978); Jastrzebska et al. (1996). However, it is by no means certain that the conductivity reported in any of these experimental studies is electronic in nature. A clear idea of the basic structural unit is critical to developing a consistent model for condensed phase charge transport in such disordered organic systems. It is also important in the context of “molecularly engineering” melanins to have the ability to create or enhance functionality in high technology applications such as bio-sensors and bio-mimetic photovoltaics Meredith and Riesz ; Riesz (2003).

It is well known that disorder plays a crucial role in determining the charge transport properties of amorphous inorganic semiconductors. In addition to the usual impurity related disorder, organic semiconductors (of which eumelanin may well be an exotic example), present several different sources of disorder associated with structural heterogeneity. For example, if one considers either the hetero-polymer or oligomer models, there are many different ways of constructing a eumelanin macromolecule based upon the monomer sequences, cross-linking positions or tautomer combinations. These possibilities have been studied theoretically for DHI-based macromolecules using a number of quantum chemical techniques. Notably, Galvao and Caldas Galvão and Caldas (1988, 1990a, 1990b) used the Hückel approximation to construct model homo-polymers. Intriguingly, they found that many of the important “semiconductor related properties of eumelanin” stabilise at a relatively small number of monomer units - five or six in fact. Although they did not consider inter-chain effects, and essentially treated chain terminations as end defects, their work was the first indication that large, extended hetero-polymeric structures are not required to explain the physical properties of melanins. Bolivar-Marinez et al. Bolívar-Marinez et al. (1999) and Bochenek and Gudowska-Nowak Bochenek and Gudowska-Nowak (2003) have used the intermediate neglect of differential overlap (INDO) and other semi-empirical methods to perform similar calculations for oligomers, whilst other authors have utilised the power of density functional theory to perform first principles calculations Stark et al. (2003); Il’ichev and Simon (2003) on single monomers or dimers.

Il’ichev and Simon Il’ichev and Simon (2003) considered the hydroquinone, (HQ) and indolequinone (IQ) forms (see figure 1), and determined them to be the most stable and hence the most likely forms of the monomers. However, several pulse radiolysis studies Lambert et al. (1989, 1990); Al-Kazwini et al. (1990, 1991); Al-Kazwini et al. (1992a, b) have indicated the presence of the semiquinone (SQ - also shown in figure 1), a tautomer of IQ, and so we also consider SQ here.

Motivated by the remarkable optical absorption characteristics of eumelanins, many of these authors have focused on simulating the absorption spectra of HQ, IQ, SQ and related small oligomers Il’ichev and Simon (2003); Bolívar-Marinez et al. (1999); Stark et al. (2003). However, knowledge of the electronic states around the highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO) levels is only part of the story. A great deal of information is also contained within the vibrational and rotational structure of the molecules. In particular, the IR and Raman spectra may be useful analytical tools. This has proven to be the case in other biomolecules van Holde et al. (1998). They can be used to confirm the accuracy of the first principles predictions by direct comparison with experimental evidence, and also to provide useful insight into local environment phenomena such as solvent-solute interactions. Furthermore, our IR and Raman spectra calculations have the advantage of coming from first principles calculations and not the semi-empirical methods often employed to calculate the absorption spectra.

In this paper we present a full density functional theory (DFT) analysis of the electronic, vibrational and rotational structure of the indolequinone (IQ), the semiquinone (SQ) and the hydroquinone (HQ). Our study is a considerable extension of previously published work since we: i) present calculated IR and Raman spectra for the first time, and ii) use the difference of self consistent fields ($`\mathrm{\Delta }`$SCF) approach to gain a more accurate understanding of the LUMO level. We find that the electronic properties of the molecules are highly dependent upon the redox form. Hence, we can say that oligomers or hetero-polymers consisting of one or more redox forms will have large variations in the local chemical potential, or, in the language of the tight binding model, strongly random site energies and possibly also random hopping integrals. This is to say that such hetero-structures will be highly disordered.

Our work is motivated by the desire to understand the implications of the redox form of the basic structural unit of eumelanins on the materials bulk properties. Such knowledge is a key starting point in any attempts to understand charge transport in and the optical properties of these disordered heterogeneous organic conductors. It also has profound implications in unravelling the mysteries of melanin biological functionality, and attempts to molecularly engineer melanin-like molecules for technological applications. These preliminary theoretical calculations are part of an ongoing quantum chemical, experimental solid-state and spectroscopic program aimed at gaining a more complete understanding of melanin structure-property-function relationships.

## II Calculation details

The chemical and electronic structures and the Raman and IR spectra were found from first principles DFT calculations. We have performed our calculations using the Naval Research Laboratory Molecular Orbital Library (NRLMOL) Pederson and Jackson (1990); Jackson and Pederson (1990); Pederson and Jackson (1991); Quong et al. (1993); Porezag and Pederson (1996); Porezag (1997); Briley et al. (1998). NRLMOL performs massively parallel electronic structure calculation using gaussian orbital methods. In particular, for a discussion of the calculation of the Raman and IR spectra see Refs. Porezag and Pederson (1996); Porezag (1997); Briley et al. (1998). Throughout we have used the Perdew, Burke and Ernzerhof (PBE) Perdew et al. (1996) exchange correlation functional, which is a generalised gradient approximation (GGA) containing no parameters. For each molecule we have fully relaxed the geometry with no symmetry constraints.

## III Results and Discussion

To benchmark our calculations we report the bond lengths and bond angles of HQ, IQ and SQ found in our calculations in table 2. Our results are in good agreement with those calculated by Bolivar-Marinez et al. Bolívar-Marinez et al. (1999) and Stark et al. Stark et al. (2003). All three molecules are planar to within numerical accuracy (c.f., Il’ichev and Simon (2003); Stark et al. (2003)).

We find, in agreement with Il’ichev and Simon Il’ichev and Simon (2003), that SQ is not the lowest energy tautomer of IQ (see table 1). Indeed we find that at 300 K there should be less than 0.1% SQ. Although we should stress that this calculation is for a single molecule *in vacuo* it is, of course, in direct contradiction to the evidence Lambert et al. (1989, 1990); Al-Kazwini et al. (1990, 1991); Al-Kazwini et al. (1992a, b) that SQ is one of the building blocks of eumelanin. However, since SQ is widely thought to play an important role in eumelanin, we also consider it here.

As well as the ground state geometries shown in figures 2, 3 and 4 we also found stable geometries of HQ and SQ with the H-O-C bond angle increased by almost 180. These alternative structures have significantly higher energies than the ground states. However, it is possible that the high energy structures may be stabilised by hydrogen bonding in a polar solvent. We therefore also considered the system $`\text{SQ}+6\text{H}_2\text{O}`$. We found that in this system the H-O-C bond is slightly larger than in SQ in vacuo (see figure 5). We also found a significant hydrogen bonded network around the H-O-C-C=O group which is dragged out of the plane of the molecule by interactions with the water molecules. This hydrogen bonding network is probably related to the fact that eumelanins are only soluble in polar solvents.

The Hückel model calculations of Galvao and Caldas Galvão and Caldas (1988) indicate that the HOMO-LUMO gap ($`\mathrm{\Delta }_{\text{HL}}`$) may be significantly larger in HQ than in either IQ or SQ. This has important implications in the context of understanding the basic eumelanin structural unit. Fundamentally, a combination of HQ and either IQ or SQ tautomers are present in an ensemble. Whether that ensemble be oligomeric or polymeric it can be expected to lead to an apparent broadening of the absorption profile. In this situation, the key question is: how many monomers of HQ, IQ and SQ are required to create the broad band UV and visible absorbance over the range 1.5 to 5 eV that is macroscopically observable for eumelanins? The answer to this question may well allow us to assess the feasibility of the various structural models. Furthermore, when one considers the possibility that molecularly engineered forms of eumelanin may be useful as functional materials in electronic devices and sensors Meredith and Riesz ; Riesz (2003), a critical design parameter is likely to be the “semiconductor” gap. It has previously been shown Galvão and Caldas (1990a) that, at least within the Hückel approximation for homopolymers, the HOMO-LUMO gap of the monomer is closely related to the HOMO-LUMO gap of the infinite polymer, which, for the polymer model in the condensed solid state, corresponds to the semiconducting gap provided there is only weak coupling between polymers. Thus controlling the HOMO-LUMO gap of macromolecules may provide a route to controlling the semiconducting gap of eumelanins. This suggestion clearly needs further investigation, but controlling the HQ to IQ or SQ ratio may be a possible route to this form of property manipulation. An important first step to investigating either of these proposals is the calculation of the HOMO-LUMO gaps of the molecules by more reliable methods than Hückel theory.

In table 3 we compare the HOMO-LUMO gap found in our calculations with that found from the Hückel method Galvão and Caldas (1988). For comparison with the Hückel calculations we need to estimate the hopping integral, $`\beta `$, in the Hückel Hamiltonian,

$`\widehat{}={\displaystyle \underset{i\sigma }{}}(\alpha +\alpha ^{}\delta _{i9})\widehat{c}_{i\sigma }^{}\widehat{c}_{i\sigma }\beta {\displaystyle \underset{ij\sigma }{}}\widehat{c}_{i\sigma }^{}\widehat{c}_{j\sigma },`$ (1)
where $`\alpha `$ is the on-site energy, $`\alpha ^{}`$ is the difference in the on-site energy between a nitrogen and a carbon atom (we take the labelling convention from figure 1 ensuring that the only nitrogen is in position $`i=9`$), $`\widehat{c}_{i\sigma }^{()}`$ annihilates (creates) an electron with spin $`\sigma `$ on site $`i`$ and $`ij`$ indicates that the sum is over nearest neighbours only. For hydrocarbons $`\beta `$ is typically of order 2.5 eV Lowe (1978). Using this value of $`\beta `$ we find that the HOMO-LUMO gaps for HQ, IQ and SQ found from the Hückel approximation are consistent with those found in our DFT calculations using the PBE functional. However, because DFT is a theory of the ground state, these calculations represent the energy gap between Kohn–Sham eigenvalues and not the true HOMO-LUMO gap of the molecules. This is known as the band gap problem Jones and Gunnarsson (1989). Additionally, it is accepted that the PBE functional can significantly underestimate the HOMO-LUMO gap. Therefore we have also employed the $`\mathrm{\Delta }`$SCF method Jones and Gunnarsson (1989) to calculate $`\mathrm{\Delta }_{\text{HL}}`$. Our results clearly reproduce the trends seen in the time dependent density functional theory (TDDFT) calculations of Il’ichev and Simon Il’ichev and Simon (2003) for IQ and SQ. The difference between the $`\mathrm{\Delta }`$SCF and TDDFT results for HQ probably indicates that the true value of $`\mathrm{\Delta }_{\text{HL}}`$ is intermediate to the values predicted by the two different methods as it is extremely unlikely that either method is incorrect by a large enough margin to allow the other to be correct.

The calculated $`\mathrm{\Delta }`$SCF HOMO-LUMO gap is also consistent with previously published semiempirical ZINDO (Zerner’s intermediate neglect of differential overlap) calculations of the optical absorption spectra Bolívar-Marinez et al. (1999); Stark et al. (2003). In particular Bolívar-Marinez et al. noted that “in the neutral state the threshold for optical absorption is around 2.0 eV for the IQ and SQ, while in the case of the HQ it is roughly 3.8 eV.” Taking these numbers as semiempirical estimates of the HOMO-LUMO gap they are in excellent agreement with our $`\mathrm{\Delta }`$SCF calculations (c.f. table 3).

It is also interesting to note that the change in the SQ conformation observed in the calculations for the $`\text{SQ}+6\text{H}_2\text{O}`$ system caused a 7% decrease in $`\mathrm{\Delta }_{\text{HL}}`$ relative to SQ in vacuo. It is possible that in a polar solvent a range of conformations exist which may lead to a broadening of the optical absorption Myers (1998).

In terms of Hückel theory, the difference between $`\mathrm{\Delta }_{\text{HL}}`$ for HQ and IQ/SQ can be understood in terms of the level of delocalisation of the HOMOs and LUMOs. It should be noted that the aromatic ring in HQ contains a resonating valence bond, whereas this feature is absent in IQ and SQ. This can be contrasted with the Hückel model result for benzene which predicts $`\mathrm{\Delta }_{\text{HL}}=2\beta `$ for the neutral molecule but if one double bond in the ring is replaced by a single bond, as in, for example, the reduced form C6H8, then $`\mathrm{\Delta }_{\text{HL}}=0`$. On this basis one might well expect the HOMO-LUMO gap of HQ to be significantly larger than those of IQ and SQ. However, detailed plots of the electron densities in the HOMO and LUMO of the three molecules (figures 2, 3 and 4) indicate that the real situation is somewhat more complicated. The electron densities are significantly more localised in the HOMOs of IQ and SQ than the electron density of the HOMO of HQ is. In particular, the electron densities in the HOMOs and LUMOs of IQ and SQ are not as similar as one might expect from figure 1. However, even in the $`\mathrm{\Delta }`$SCF calculation, the HOMO-LUMO gap is a factor of two larger in HQ than it is in SQ and IQ.

In the absence of knowledge of the chemical structure of the eumelanins it is difficult to predict how a change in the HOMO-LUMO gap of the monomer will effect the electronic structure. However it is likely that changing $`\mathrm{\Delta }_{\text{HL}}`$ by a factor of two will have a dramatic effect on the electronic structure of the melanin macromolecule. For example, it is reasonable to expect that the HOMO-LUMO gap of a macromolecule will be related to $`\mathrm{\Delta }_{\text{HL}}`$ (c.f. Ref. Galvão and Caldas (1990a)). Thus we expect that the HOMO-LUMO gap for a macromolecule of HQ will be significantly larger than that of a macromolecule of IQ or SQ. For a macromolecule containing both HQ and, say, IQ this change in the electronic structure will act as a source of disorder and thus have a dramatic effect on the electronic transport properties.

It is also difficult, without detailed knowledge of the chemical structure of the eumelanins, to predict exactly how the variation of $`\mathrm{\Delta }_{\text{HL}}`$ affects the optical absorption. However, it seems possible that the range of $`\mathrm{\Delta }_{\text{HL}}`$ between different molecules may produce the continuum of HOMO-LUMO gaps in eumelanin macromolecules and thus play an important role in explaining the observed the broad band optical absorption.

It is therefore important to be able to identify the monomeric content of a sample of eumelanin, ideally this should be done in situ and non-destructively. Motivated by this fact we have calculated both the Raman and infrared (IR) spectra of IQ, SQ and HQ, (see tables 4, 5 and 6). This has several uses. Firstly, the prediction allows for experimental testing of the accuracy of our calculations. Secondly, it can be seen that there are notable differences in both the Raman and IR spectra of the three monomers (see, for example, figures 6 and 7). The Raman and IR spectra could therefore be used for the in situ, non-destructive identification of the monomeric content of macromolecules. This could be particularly valuable in the engineering of devices from eumelanins, as, by varying the ratio of indole-quinones, it may be possible to achieve control of the band-gap of the material. This may also be a useful analytical tool for the analysis of the structure of both synthetic and naturally occurring melanin, and thus be helpful in determining the structure-property-function relationships that control their behaviour.

Comparison of our calculated IR and Raman spectra with standard tables for organic molecules Nakanishi and Solomon (1977); Schrader (1995) shows that the calculated phonon spectra are in broad agreement with known IR absorption band values (see table 7). It can also be seen from table 7 that the differences in the spectra, specifically which vibrations are absent for a given molecule can be understood in terms of their structure, c.f. figure 1, as expected.

It is important to stress that our calculated IR and Raman spectra are for gaseous phase monomers. Thus, the comparison of such results with experiment may be complicated by several factors. Firstly the formation of macromolecules may result in a significant change in the spectra. However, it is known that in other biomolecules the effect of polymerisation only slightly broadens the IR and Raman spectra, for example the spectra of simple amino acids is only slightly broadened in proteins Xie et al. (2000). Secondly solvent effects may largely ‘wash out’ the individual features of the spectra Reichardt (1998); Myers (1998). It is clear from our study of the $`\text{SQ}+6\text{H}_2\text{O}`$ system that there will be significant hydrogen bonding in the solvent, hydrogen bonding can shift the phonon frequencies van Holde et al. (1998). Therefore it remains to be seen if, experimentally, the resolution is good enough to resolve the differences in the spectra of the monomers. However, an obvious way to avoid solvent effects is to conduct the experiments in the solid state.

## IV Conclusions

We have carried out first principles density functional calculations for the hydroquinone (HQ), the indolequinone (IQ) and the semiquinone (SQ). The calculated gaseous phase structure is in good agreement with previous calculations. We have used the $`\mathrm{\Delta }`$SCF method to study the HOMO-LUMO gap. Our results are consistent with the threshold for optical absorption found by semiempirical methods. Specifically, we found that the HOMO-LUMO gap is similar in IQ and SQ but approximately twice as large in HQ. The possibility of using this difference in the HOMO-LUMO gap to engineer the electronic properties of eumelanins at the molecular level has been discussed. We have also suggested that the difference in the HOMO-LUMO gap of the different monomers could lead to a large range of HOMO-LUMO gaps in eumelanin macromolecules and thus be related to the observed broad band optical absorption.

As the structure of macromolecules of these monomers is not known we have also calculated the IR and Raman spectra of the three monomers from first principles. A comparison of these results with experiment would represent a stringent test of the density functional calculations as these calculations do not rely on additional semiempirical information. It was shown that the IR and Raman spectra have potential as a analytical tool for these materials. Each of the monomers have significantly different spectral signatures, therefore the IR or Raman spectra could be used in situ to non-destructively investigate the monomeric content of macromolecules. It is hoped that this may be helpful in determining the structure of both natural and synthetic eumelanins, and hence aid our attempts to understand their biological functionality.

## Acknowledgements

We would like to thank John Dobson, Joel Gilmore, Peter Innis, Urban Lundin, Jeffrey Reimers, Jennifer Riesz Cathy Stampfl and Linh Tran for helpful discussions. The work at the University of Queensland was supported by the Australian Research Council. BJP was supported in part by US Navy Grant N00014-03-1-4115. KB was supported in part by US Navy Grant N00014-03-1-4116. TB was supported in part by US Navy Grant N00014-02-1-1046. MRP was supported in part by ONR and the DoD HPC CHSSI initiative. BJP, RHM and MRP appreciate the hospitality of the Kavli Institute of Theoretical Physics, where this work was initiated. The KITP is supported by NSF through grant PHY99-07949.


text:
# Thermodynamic and transport properties of underdoped cuprates from ARPES data

## 1 Introduction

Thermodynamics and the transport properties of high-$`T_\text{c}`$ cuprates are highly affected by the strong electron correlation near the filling-control Mott transition. One of the most unconventional thermodynamic properties is that the electronic specific heat coefficient $`\gamma `$ shows pseudogap phenomena in the underdoped region . The electronic resistivity show unconventional features such as $`T`$-liner temperature dependence in the optimally doped region and metallic behavior well exceeding the Ioffe-Regel limit in the lightly doped region . In order to elucidate the origin of these unconventional phenomena, angle-resolved photoemission spectroscopy (ARPES) is an extremely powerful technique because one can directly observe the spectral function $`A(𝐤,\omega )`$, i.e., the electronic structure, near the Fermi level ($`E_\mathrm{F}`$). Therefore, information obtained from ARPES can be directly compared with thermodynamic properties such as the electronic specific heat. Also, the ARPES spectra near $`E_\mathrm{F}`$ give significant information about the transport properties. In this article, we shall attempt to understand the unconventional thermodynamics and transport properties based on ARPES spectra mostly just at $`E_\mathrm{F}`$ .

## 2 Thermodynamic properties

The electronic specific heat coefficient $`\gamma `$ of a Fermi liquid is given by $`\gamma =\pi ^2k_B^2/3N(0)^{}`$ , where $`N(0)^{}`$ is the quasi-particle (QP) density at $`E_\mathrm{F}`$. In a single-band system like single-layer cuprates, each momentum contributes one QP, if the system is a Fermi liquid. Therefore, if one successfully fits the dispersion of the QP band to some energy band model such as the tight-binding (TB) model, the density of states (DOS) of that band yields $`N(0)^{}`$.

In Fig. Thermodynamic and transport properties of underdoped cuprates from ARPES data(a), the DOS determined by a TB fit to the ARPES results of La2-xSrxCuO4 (LSCO) is compared with the specific heat coefficient $`\gamma `$ of LSCO samples where the $`T_\text{c}`$ is suppressed by Zn-doping . The TB parameters have been determined to reproduce the experimentally observed Fermi surface and the slope of the band dispersions in the nodal (0,0)-($`\pi `$,$`\pi `$) direction. Here, fine structures of “kink” in the energy dispersions  are neglected in the TB fit. The TB band DOS has a maximum around $`x`$=0.15 because the flat band near ($`\pi ,0`$) becomes closest to the Fermi level. If the system is a Fermi liquid, the TB band DOS and $`\gamma `$ should be the same. However, the $`\gamma `$ decreases much faster than TB band DOS with decreasing $`x`$. This difference indicates that the QP density at $`E_\mathrm{F}`$ that would be expected for a normal Fermi liquid is indeed depleted in the underdoped materials and demonstrates how the opening of the pseudo-gap occurs around ($`\pi ,0`$).

The doping dependence of the DOS of angle integrated photoemission (AIPES) spectra at $`E_\mathrm{F}`$ are also shown in Fig. Thermodynamic and transport properties of underdoped cuprates from ARPES data(a). AIPES spectra have been obtained by integrating the ARPES spectra with respect to $`𝐤`$ over the Brillouin zone. As a whole, the AIPES DOS decreases with decreasing $`x`$ for $`x<0.22`$ due to the pseudo-gap opening. However, a closer inspection reveals that the AIPES DOS decreases faster than $`\gamma `$ with decreasing $`x`$. Note that the DOS at $`E_\mathrm{F}`$ measured by AIPES is given by $`zN(0)^{}`$, where $`z`$ is the renormalization factor. The faster decrease of the AIPES intensity with decreasing $`x`$ is due to the effects of $`z`$, which also decreases with decreasing $`x`$ . From the above comparison, one can say that the AIPES DOS reflects the effects of both the pseudogap opening and the spectral weight renormalization $`z`$.

Relationship between the ARPES spectral weight and the superfluid density would also give important key about the superconductivity . The doping dependence of the ARPES spectral weight at $`E_\mathrm{F}`$ is compared with that of the penetration depth $`1/\lambda ^2\rho _s/m^{}`$  in Fig. Thermodynamic and transport properties of underdoped cuprates from ARPES data(b). The spectral weight at angle $`\alpha `$ (see inset) is defined by the integrated spectral weight of the spectral function $`A(𝐤,\omega )`$ at $`E_\mathrm{F}`$ over the momentum.

The decrease of both the in-plane and out-of plane superfluid densities for $`x<`$0.20 is due to the pseudogap opening near ($`\pi `$,0). Particularly, the doping dependence of the out-of-plane $`(\rho _s/m^{})_c1/\lambda _c^2`$ well corresponds to that of the spectral weight around $``$($`\pi `$,0) ($`\alpha `$=0). This is because the out-of-plane $`(\rho _s/m^{})_c1/\lambda _c^2`$ has large matrix element around ($`\pi `$,0). On the other hand, the nodal direction is expected to be the main contribution to in-plane $`(\rho _s/m^{})_{ab}1/\lambda _{ab}^2`$ because of the lighter effective mass. However, the decrease with decreasing $`x`$ is somewhat different from the doping dependence of the nodal spectral weight ($`\alpha `$=45). This implies that not only the nodal QP’s but also QP’s over an extended $`k`$-region of the “Fermi arc” condense into the superfluid.

## 3 Transport properties

The transport properties of high-$`T_\text{c}`$ cuprates show unconventional behavior reflecting the strongly electron correlation. Because of the pseudo-gap opening around ($`\pi `$,0), the electronic states around the node mainly contribute to the in-plane transport properties in the underdoped region. This portion of the Fermi surface survives in the upderdoped compounds as a “Fermi arc” . We shall estimate the in-plane electrical resistivity from the observed “Fermi arc” using the conventional Drude formula $`\rho =m^{}/ne^2\tau `$. Here, the renormalized mass $`m^{}=\mathrm{}k_\mathrm{F}/v_\mathrm{F}`$ can be determined from the quantities measured by ARPES, namely, the Fermi velocity $`v_\mathrm{F}`$ and the Fermi momentum $`k_\mathrm{F}`$ of the hole-like Fermi surface centered at ($`\pi ,\pi `$). From the MDC width $`\mathrm{\Delta }k`$ at $`E_\mathrm{F}`$, the mean free path is obtained as $`l=1/\mathrm{\Delta }k`$, yielding the mean free scattering time of $`\tau =l/v_\mathrm{F}`$. Thus, one can evaluate $`\rho =m^{}/ne^2\tau =\mathrm{}k_\mathrm{F}\mathrm{\Delta }k/ne^2`$ from ARPES spectra.

Figure Thermodynamic and transport properties of underdoped cuprates from ARPES data(a) shows an ARPES intensity plot in the $`E`$-$`k`$ space for an LSCO ($`x`$ = 0.03) sample in the (0,0)-($`\pi ,\pi `$) direction. The energy dispersion determined by the peaks in the MDC shows a “kink” feature around 70 meV below $`E_\mathrm{F}`$, which is presumably due to electron-phonon coupling . $`v_\mathrm{F}1/m^{}`$ in the node direction is almost doping independent $`v_\mathrm{F}1.8`$ eV Å for LSCO .

Now, we discuss how to estimate the carrier number $`n`$. Although the Fermi-liquid picture yields the Fermi surface volume $`1x`$, Hall coefficient results show unconventional behavior $`nx`$ in the underdoped region . The observed “Fermi arc” picture gives two possibilities to explain the evolution of the carrier number like $`nx`$ . One possibility is that the length of the “Fermi arc” increases with $`x`$ and yields $`nx`$. However, the “arc” length estimated from the ARPES intensity does not change enough to explain the $`nx`$ behavior. Another possibility is the evolution of the nodal spectral weight with $`x`$. Indeed, the spectral weight as a function of $`x`$ is remarkably similar to the $`x`$-dependence of the Hall coefficient .

Figure Thermodynamic and transport properties of underdoped cuprates from ARPES datab shows the doping dependence of the estimated in-plane electrical resistivity $`\rho _{\mathrm{ARPES}}`$ based on the Drude formula by using the parameters $`k_\mathrm{F}`$, $`\mathrm{\Delta }k`$ and $`n`$ as described above. ARPES intensity $`n_{\mathrm{PES}}n_{\mathrm{Hall}}`$, $`x`$ (“small Fermi surface volume”) and $`1x`$ (“large Fermi surface volume”) are considered as the carrier number $`n`$. $`\rho _{\mathrm{ARPES}}`$ evaluated using $`n=x`$ and $`n_{\mathrm{PES}}`$ are in good agreement with the transport results in the underdoped region. On the other hand, $`\rho _{\mathrm{ARPES}}`$ evaluated using $`n=1x`$ is by more than one order of magnitude smaller than the transport data in the underdoped region. This suggests that the pseudogap and the small renormalized spectral weight is the origin of the small carrier number and the high electrical resistivity for the underdoped samples. On the other hand, the discrepancy between $`\rho _{\mathrm{ARPES}}`$ ($`n=x`$,$`n_{\mathrm{PES}}`$) and the transport data becomes pronounced in the optimally doped region. This discrepancy may indicate a cross-over from the Fermi arc to the large Fermi surface and/or may be due to the contribution from the ($`\pi `$,0) region, which is not precisely taken into account in the present simplified analysis. We leave such quantitative analysis for future investigation. As a whole, the estimated $`\rho _{\mathrm{ARPES}}`$ for $`n_{\mathrm{PES}}`$ explains the doping dependence of the transport data of the underdoped LSCO qualitatively well.

We also show the case of YBa2Cu3Oy(YBCO) in Fig. Thermodynamic and transport properties of underdoped cuprates from ARPES data(a)(b). Relationship between the oxygen content $`y`$ and the hole concentration $`\delta `$ has been estimated from the electrical resistivity and the thermopower . The MDC peak width of YBCO is about twice as large as that of LSCO and the MCD line shape of YBCO could not be well fitted to a single Lorentzian. We consider that this is due to the bilayer splitting in the nodal direction as predicted by band-structure calculations . We have therefore fitted the MDC to two Lorentzians to estimate $`\mathrm{\Delta }k`$, assuming that the chain band does not exist at $`E_\mathrm{F}`$ for such low doping. Fig.Thermodynamic and transport properties of underdoped cuprates from ARPES data(b) shows the estimated value of $`\rho `$ in the same way as Fig.Thermodynamic and transport properties of underdoped cuprates from ARPES data(b). While the transport data  is closer to $`\rho _{\mathrm{ARPES}}(n=\delta )`$ for $`\delta =0.04`$, it approaches $`\rho _{\mathrm{ARPES}}(n=1\delta )`$ with increasing $`\delta `$, again implying the crossover from the Fermi arc to the large Fermi surface.

Because the Fermi velocity $`v_\mathrm{F}`$ does not seem to change with temperature, the temperature dependence of the inverse life time $`1/\tau `$ is determined by that of $`\mathrm{\Delta }k=1/l=1/v_\mathrm{F}\tau `$. Figure Thermodynamic and transport properties of underdoped cuprates from ARPES data(a) shows the temperature dependence of MDC at $`E_\mathrm{F}`$ in the nodal direction for the $`x`$=0.07 sample. The MDC width $`\mathrm{\Delta }k`$ is systematically broadened with increasing temperature, indicating a decrease of the mean-free path. As shown in Fig.Thermodynamic and transport properties of underdoped cuprates from ARPES data(b), the temperature dependence of $`\mathrm{\Delta }k`$ well scales with that of the electric resistivity, indicating that the temperature dependence of $`\mathrm{\Delta }k`$ determines that of the resistivity. In optimally doped Bi2212, the linear temperature dependence of the MDC width which scales with resistivity was observed and discussed in the context of marginal Fermi-liquid picture . The present result also shows a $`T`$-linear behavior in $`\mathrm{\Delta }k`$ over a wide temperature range. However, considering that electron-phonon interaction causes the “kink” behavior  and contribute to the low energy properties, the MDC width would be affected by phonons.

## 4 Conclusion

In summary, we have discussed the relationship between the ARPES results for high-$`T_\text{c}`$ cuprates and their thermodynamics and transport properties. From the comparison of the ARPES and $`\gamma `$ and $`\rho _s`$, the pseudogap feature and the spectral weight renormalization are discussed. We have estimated the electric resistivity from the observed spectra and have demonstrated that the obtained results explain the transport properties. The temperature dependence of the MDC width $`\mathrm{\Delta }k`$ in the nodal direction well scales with that of the electric resistivity.

We are grateful to H. Fukuyama for enlightening discussions. This work was supported by a Grant-in-Aid for Scientific Research “Novel Quantum Phenomena in Transition Metal Oxides” from the Ministry of Education, Science, Culture, Sports and Technology of Japan and the New Energy and Industrial Technology Development Organization (NEDO). Experimental data were recorded at ALS and SSRL which are operated by the Department of Energy’s Office of Basic Energy Science, Division of Material Sciences. SSRL is also operated by Division of Chemical Sciences.

Fig. 1 (a) Comparison the electronic specific heat coefficient $`\gamma `$  and the TB band DOS at $`E_F`$. The doping dependence of AIPES DOS at $`E_\mathrm{F}`$ is also shown. (b) Doping dependence of the spectral weight around nodal and ($`\pi `$,0) regions. The in-plane and out-of-plane superfluid densities deduced from the magnetic penetration depth ($`\lambda _{ab}`$, $`\lambda _c`$)  are compared with them.

Fig. 2 (a) Spectral intensity in the $`E`$-$`k`$ space along the nodal cut and the MDC at the Fermi level. (b) Comparison of the electrical resistivity at 20 K  and the resistivity estimated from ARPES. The carrier number $`n`$ is assumed to be $`n=x,n=1x`$ and $`n=n_{\mathrm{PES}}`$.

Fig. 3 (a) Spectral intensity of YBCO in the $`E`$-$`k`$ space along the nodal cut and the MDC at the Fermi level. (b) Comparison of the electrical resistivity at 10 K  and the resistivity estimated from ARPES. The carrier number $`n`$ is assumed to be $`n=\delta `$ and $`n=1\delta `$.

Fig. 4 (a) Temperature dependence of the MDC at $`E_\mathrm{F}`$ in the nodal direction for $`x`$=0.07. (b) Comparison of the temperature dependence of the MDC momentum width $`\mathrm{\Delta }k`$ with that of the electrical resistivity.

Fig. 1

Fig. 2

Fig. 3

Fig. 4


text:
# Extended moment formation and magnetic ordering in the trigonal chain compound 𝐂𝐚_𝟑𝐂𝐨_𝟐𝐎_𝟔

## 1 Introduction

Low-dimensional physical systems are attracting a lot of attention since long due to the occurence of exciting physical properties deviating distinctly from those known from three-dimensional systems. This includes both, new types of ordering phenomena being due to, e.g., Fermi-surface instabilities and new types of excitations as resulting, e.g., from spin-charge separation. Furthermore, being embedded in three-dimensional crystals, low-dimensional atomic arrangements give rise to different length scales and coupling strengths and thus may lead to rich phase diagrams. Finally, the existence of magnetic moments in such systems raises questions concerning the relevant exchange mechanisms and their anisotropies.

In recent years, interest especially in one-dimensional systems has focused on a new class of transition metal oxides of the general formula $`\mathrm{A}_3^{}\mathrm{ABO}_6`$ ($`\mathrm{A}^{}=\mathrm{Ca},\mathrm{Sr},\mathrm{Ba}`$; $`\mathrm{A},\mathrm{B}=`$ transition metal), which crystallize in the trigonal $`\mathrm{K}_4\mathrm{CdCl}_6`$ structure . The space group is $`R\overline{3}c`$ ($`D_{3d}^6`$, No. 167). In these compounds, transition metal-oxygen polyhedra form well separated chains running parallel to the trigonal axis; see, e.g., Refs.  for a representation of the crystal structure. Space between the chains is filled with the $`\mathrm{A}^{}`$ cations. Each chain consists of alternating, face-sharing $`\mathrm{AO}_6`$ trigonal prisms and $`\mathrm{BO}_6`$ octahedra. Typically, the ratio of interchain to intrachain metal-metal distance is of the order of two, this fact explaining the pronounced one-dimensionality. The chains themselves are arranged on a triangular lattice. As a consequence, in addition to showing the abovementioned unique properties these compounds allow to study geometric frustration effects as partial disorder and spin-glass like behaviour.

Prominent members of this class are the cobaltates $`\mathrm{Ca}_3\mathrm{CoBO}_6`$ with $`\mathrm{B}=\mathrm{Co}`$, $`\mathrm{Rh}`$, $`\mathrm{Ir}`$, which have gained much interest due to their exciting magnetic properties. Here we concentrate on $`\mathrm{Ca}_3\mathrm{Co}_2\mathrm{O}_6`$, which has two inequivalent cobalt sites, one in octahedral environment (labelled Co1) and the other (Co2) centered in the trigonal prisms . According to powder neutron diffraction data reported by Fjellvåg et al. room-temperature lattice constants amount to $`a_{hex}=9.0793`$ Å and $`c_{hex}=10.381`$ Å. Data taken at 10 K by Aasland et al. resulted in $`a_{hex}=9.060`$ Å and $`c_{hex}=10.366`$ Å with only slight changes of the atomic positional parameters as compared to 298 K . The latter data were used in the calculations.

While the Co–Co distance within the chains amounts to 2.595 Å, next-nearest neighbour distances across the chains are 5.31 Å . This value is somewhat larger than the interchain separation of 5.24 Å, since neighbouring chains are shifted parallel to the trigonal axis by 1/3 of the intrachain separation. According to the neutron diffraction data, the octahedral Co1-O distances of 1.916 Å are considerably shorter than the Co2-O distances, which amount to 2.062 Å. From the difference it was expected that $`\mathrm{Ca}_3\mathrm{Co}_2\mathrm{O}_6`$ might show charge or spin order at low temperatures . Finally, the short Co1–Co2 distance within the chains, which is a result of the face-sharing of the polyhedra and is close to the value of 2.51 Å in metallic cobalt, has been taken as indicative for strong metal-metal bonding .

Magnetic susceptibility data revealed the onset of magnetic ordering below $`\mathrm{T}_{\mathrm{C1}}=24`$ K . In addition, a second transition was found at $`\mathrm{T}_{\mathrm{C2}}=12`$ K . From Curie-Weiss behaviour above 80 K an effective magnetic moment of $`5.7\mu _B`$ and a paramagnetic Curie temperature of 28 K pointing to predominant ferromagnetic ordering was inferred. Neutron diffraction measurements by Aasland et al. revealed local magnetic moments of $`0.08\mu _B`$ and $`3.00\mu _B`$ for octahedral Co1 and trigonal prismatic Co2, respectively . All moments are aligned along the trigonal axis. From these data it was furthermore concluded that the coupling is ferromagnetic within the chains and antiferromagnetic across the chains . While $`\mathrm{T}_{\mathrm{C1}}`$ has been associated with the ferromagnetic intrachain coupling, $`\mathrm{T}_{\mathrm{C2}}`$ would correspond to long-range three-dimensional ordering .

Magnetization measurements showed a plateau at $`1.31\mu _B`$ per f.u. for low field and a steep increase by a factor of three at about 3.5 T, which was interpreted as a ferri- to ferromagnetic transition . However, the final magnetic moment of about $`4\mu _B`$ was difficult to reconcile with the values found for the local cobalt moments . While Aasland et al. proposed a magnetic structure with 1/3 of the chains being antiparallel to the remaining 2/3 , Kageyama et al. argued that at low field and temperatures between 10 and 24 K the system assumes a partially disordered antiferromagnetic (PDA) state . The existence of this state in different materials is still a matter of controversal discussion. While it was also reported for $`\mathrm{Ca}_3\mathrm{CoRhO}_6`$ , recent magnetization and heat capacity data allowed to exclude the PDA state for $`\mathrm{Ca}_3\mathrm{CoIrO}_6`$ . For $`\mathrm{Ca}_3\mathrm{Co}_2\mathrm{O}_6`$ it was shown that pressure leads to increase of $`\mathrm{T}_{\mathrm{C2}}`$ and, hence, stabilization of the ferrimagnetic phase . Magnetization data taken below $`\mathrm{T}_{\mathrm{C2}}`$ on single crystals by Maignan et al. displayed a sequence of up to four different magnetization plateaus pointing to a variety of competing magnetic structures . Moreover, AC-$`\chi `$ measurements by Maignan et al. showed a remarkably strong frequency shift, which has been interpreted as a measure of the frustration between different magnetic configurations . These results supported interpretations of $`\mathrm{Ca}_3\mathrm{Co}_2\mathrm{O}_6`$ as a triangular lattice of antiferromagnetically coupled Ising spins, each formed from a single chain. However, the exact nature of the magnetic ground state of this material is still unclear.

In this letter we report on electronic structure calculations for $`\mathrm{Ca}_3\mathrm{Co}_2\mathrm{O}_6`$, which include the spin-degenerate case as well as spin-polarized configurations assuming both the spin-structure proposed by Aasland et al. and the hypothetical case of ferromagnetic alignment of the chains. Our calculations reveal i) low- and high-spin moments, respectively, at the octahedral and trigonal prismatic cobalt sites, ii) a rather large oxygen moment due to polarization by the high-spin cobalt sites, which together give rise to the formation of extended but still well localized $`\mathrm{CoO}_6`$ moments and resolve the discrepancy between the magnetization and neutron diffraction data, and iii) coupling of these extended moments by ferromagnetic exchange via the $`3d`$ states of the low-spin octahedral cobalt atoms.

## 2 Methodology

The calculations were performed using the scalar-relativistic augmented spherical wave (ASW) method . In order to represent the correct shape of the crystal potential in the large voids of the open crystal structure, additional augmentation spheres were inserted. Optimal augmentation sphere positions as well as radii of all spheres were automatically generated by the sphere geometry optimization (SGO) algorithm . Self-consistency was achieved by an efficient algorithm for convergence acceleration . Brillouin zone sampling was done using an increased number of $`𝐤`$-points ranging from 28 to 408 points within the irreducible wedge.

## 3 Results and Discussion

### 3.1 Spin-Degenerate Calculations

In a first step the electronic properties of $`\mathrm{Ca}_3\mathrm{Co}_2\mathrm{O}_6`$ were calculated with spin-degeneracy enforced. The resulting electronic structure and partial densities of states (DOS) are displayed in Fig. 1.

All bands display rather strong dispersion parallel to the chain axis, i.e. along the line $`\mathrm{\Gamma }`$-A, reflecting the one-dimensionality of the compound. However, the finite perpendicular dispersion points to albeit weak interchain coupling.

Five groups of bands are identified. In the energy range from -6.4 to -2.0 eV 36 bands are observed, which trace back mainly to the O $`2p`$ states. In contrast, the next three groups of bands, which extend from -1.7 eV to the valence band maximum, from 0.6 to 1.0 eV, and from 1.8 to 2.2 eV, comprise twelve and two times four bands, respectively, and derive mainly from the Co $`3d`$ states. Finally, bands starting above 4 eV are of Co $`4s`$ character. As is obvious from the partial DOS, $`p`$$`d`$ hybridization causes substantial $`d`$ and $`p`$ contributions, respectively, below and above -2 eV reaching about 40% especially in the upper valence band and the conduction bands. Whereas $`d`$ contributions from the octahedrally coordinated Co1 show up mainly below -3.8 eV and give rise to the peak at -4.8 eV, the Co2 $`d`$ contributions are equally spread over the whole energy range between -6 and -2 eV. In passing we mention the $`d`$ admixture to the peaks at 0.8 and 2.0 eV, which is of almost pure Co2 and Co1 character, respectively. From the partial DOS the cobalt atoms can be formally assigned a $`d^6`$ configuration, which, together with the absence of charge disproportionation, is in perfect agreement with the neutron scattering results .

Crystal field splitting due to the octahedral and trigonal prismatic environment of Co1 and Co2 atoms, respectively, is revealed on closer inspection of the partial $`d`$ DOS as displayed in Fig. 2.

Note that we used for both sites the same global coordinate system with the $`z`$ axis parallel to the chain axis. It deviates from the natural coordinate system for octahedral coordination with the Cartesian axes parallel to the metal-oxygen bonds due to rotation of the octahedra centered at the Co1 sites. As a consequence, both the $`t_{2g}`$ and $`e_g`$ bands arise as admixtures of all five $`d`$ states (referred to the global coordinate system), a situation well known from other trigonal and hexagonal systems as, e.g., $`\mathrm{V}_2\mathrm{O}_3`$ . Nevertheless, at the Co1 sites, the octahedral crystal field leads to nearly perfect splitting into fully occupied $`t_{2g}`$ and empty $`e_g`$ states. In contrast, the trigonal crystal field at the Co2 sites leads to splitting into non-degenerate $`d_{3z^2r^2}`$ as well as doubly degenerate $`d_{xy,x^2y^2}`$ and $`d_{xz,yz}`$ states. While the latter appear almost exclusively in the peak around 0.75 eV, the former three orbitals dominate the upper valence band.

Apart from the abovementioned strong $`p`$-$`d`$ hybridization two findings are important for the understanding of the magnetic properties. First, as has already been mentioned, the lower and upper conduction band, respectively, centered at 0.8 and 2 eV are due exclusively to the Co2 $`d_{xz,yz}`$ and Co1 $`e_g`$ states. This is related to the fact that the latter states experience $`\sigma `$-type $`p`$-$`d`$ bonding and, being antibonding, are thus pushed to higher energies whereas the $`p`$-$`d`$-bonding of the $`d_{xz,yz}`$ orbitals is less strong. In the present situation this difference lets us expect that in a spin-polarized calculation the Co2 sites will carry larger magnetic moments. Second, we point to the large width of the $`d_{3z^2r^2}`$ bands as obvious from both the Co1 and Co2 partial DOS, which is indicative of the strong metal-metal overlap along the chain axis. This has been confirmed by the calculated crystal orbital overlap population (COOP) and covalent bond energies as well as a more detailed analysis of the electronic wave functions. According to the latter the $`d_{3z^2r^2}`$ states give rise to the dispersionless bands along $`\mathrm{\Gamma }`$-M-K-$`\mathrm{\Gamma }`$ at -1.5 eV and at $`\mathrm{E}_\mathrm{F}`$ and show finite dispersion only parallel to $`\mathrm{\Gamma }`$-A. Within these bands the bonding states are more Co2-like, whereas the antibonding states have increased Co1-like character.

### 3.2 Spin-Polarized Calculations assuming Ferromagnetic Order

In a second step, spin-polarized calculations were performed. Yet, in view of the predominant influence of the ferromagnetic interchain coupling as deduced from the experimental data we started out from the hypothetical situation, where all chains are aligned in parallel, before turning to the ferrimagnetic spin-structure proposed by Aasland et al. . Insofar our procedure is along the same line of reasoning as the calculations by Whangbo et al., who, however, did not deal with the observed ferrimagnetic structure and concentrated mainly on the origin of the cobalt moments .

For the assumed ferromagnetic order we find a converged solution with well localized magnetic moments of 2.73, 0.35, and 0.14 $`\mu _B`$, respectively, on Co2, Co1, and oxygen. These values reflect the experimental finding of a high and low spin state of the cobalt atoms at the trigonal prismatic and octahedral sites, respectively, and are close to those reported by Whangbo et al. as well as the neutron diffraction results of $`3.00`$ and $`0.08\mu _B`$ for the cobalt atoms. The total moment per unit cell, i.e. per two formula units, amounts to exactly $`8.0\mu _B`$. Worth mentioning are the rather high magnetic moments associated with the oxygen sites. Summing up to almost 1 $`\mu _B`$ per formula unit, they provide an easy explanation for the discrepancy between the neutron scattering data and the magnetization in high field.

The Co $`3d`$ partial DOS given in Fig. 3

display a rather complex relationship between the electronic and magnetic properties. Due to the reduced magnetic moment the Co1 $`3d`$ states show only small spin splitting of the order of 0.3 eV, hence, much smaller than the crystal field splitting. As a consequence, the magnetic moment is carried alone by the $`t_{2g}`$ states. Within this group of bands spin-majority $`d_{3z^2r^2}`$ states experience a pronounced localization, whereas the spin-minority states have a much larger band width. For the remaining $`t_{2g}`$ states things are reversed. Co1 $`d`$ contributions below -3.5 eV display small spin-splitting but give only a negligible contribution to the magnetic moment. This is related to the fact that $`d`$-$`p`$ overlap within the octahedra is strongest for states of $`e_g`$ symmetry, which are completely occupied or empty. The situation is different at the trigonal prismatic sites, where the spin-splitting is of the same order of magnitude as the crystal field splittings. As a consequence, spin-majority states are shifted to the energy region, where the oxygen $`2p`$ states dominate. At the same time, these states experience considerable delocalization. In contrast, spin-minority states form sharp peaks at and above the Fermi energy but give only minor contributions to the occupied states.

The different distributions of the electronic states in the octahedra and trigonal prisms have interesting consequences for the magnetic moments. This becomes clearer from Fig. 4,

which displays energy and orbital-dependent spin magnetizations, i.e. the integrated differences of the spin-dependent partial DOS. For the Co1 sites it is easily seen that the small magnetic moment is mainly carried by the $`d_{3z^2r^2}`$ states. It starts to build up in the energy region above -1.5 eV and is thus due to orbitals well localized within the Co1 atoms. This is completely different for the trigonal prismatic sites. There the magnetic moment is rather isotropic, i.e. it grows out of similar amounts from all five $`d`$ states. According to Fig. 4 the $`d_{3z^2r^2}`$ orbitals contribute $`0.6\mu _B`$, while the doubly degenerate $`d_{xz,yz}`$ and $`d_{x^2y^2,xy}`$ states give $`1.3\mu _B`$ and $`0.8\mu _B`$, respectively. Furthermore, as is clearly seen in Fig. 4, all these moments start to build up already well within the energy region, where the oxygen states are dominating, thus reflecting the considerable oxygen polarization. To conclude, the strong $`p`$-$`d`$ hybridization within the trigonal prisms, which has a strength between the rather nonbonding $`\pi `$-type and the directed $`\sigma `$-type overlap of the octahedra, leads to the formation of a combined isotropic Co2 $`d`$-O $`p`$ magnetic moment. Note that the oxygen $`p`$ orbitals participate in the local moments and not in the exchange interaction between the prisms. This situation goes under the name extended moment formation and has been previously observed in copper oxides . Finally, coupling between these extended localized moments is mediated by ferromagnetic exchange interaction via the $`d_{3z^2r^2}`$ orbitals of the low-spin Co1 atoms.

### 3.3 Spin-Polarized Calculations with Ferrimagnetic Order

In a last step, we turn to the ferrimagnetic structure, where 1/3 of the chains have their spins antiparallel to the remaining chains. The results obtained from these calculations are very similar to those for the hypothetical ferromagnetic structure and, hence, provide an a posteriori justification for the validity of the latter. Local magnetic moments of 2.72, 0.34, and 0.14 $`\mu _B`$, respectively, are found for Co2, Co1, and oxygen and the total moment per unit cell again amounts to $`8.00\mu _B`$. According to the calculated magnetic moments as well as the partial DOS, which likewise look very similar to those obtained for the ferromagnetic situation, the differences between the spin up and down chains are rather small. Yet, slightly sharper peaks, i.e. stronger localization of the Co $`3d`$ states, are observed for the minority chains. The general situation obtained for the ferrimagnetic structure is thus very similar to that grown out of the calculations for ferromagnetically coupled chains.

## 4 Conclusion

In conclusion, electronic structure calculations for the magnetic chain compound $`\mathrm{Ca}_3\mathrm{Co}_2\mathrm{O}_6`$ have revealed strong influence of the local coordination polyhedra on the electronic and magnetic properties. In particular, shorter bond lengths and strong $`\sigma `$-type bonding in the octahedra allow for only small spin-splitting at these sites. In contrast, due to smaller crystal field splittings in trigonal prismatic surrounding the Co2 atoms are found in a high-spin state. Strong $`d`$-$`p`$ hybridization causes rather large oxygen magnetic moments, which, together with the cobalt moments at the trigonal prismatic sites, take part in the formation of extended but still well localized magnetic moments. Moreover, intrachain exchange coupling between these effective moments is based on metal-metal overlap of the cobalt $`d_{3z^2r^2}`$ orbitals along the chains and mediated by ferromagnetic exchange interaction via the $`d_{3z^2r^2}`$ orbitals of the low-spin cobalt atoms. Recently, mapping these results onto a Heisenberg model we were able to underline the importance of the interplay of different crystal field splittings, $`d`$-$`p`$ hybridizations, and metal-metal overlap for the ferromagnetic intrachain order . Finally, according to a detailed analysis of the oxygen magnetic moments as well as the analysis within the Heisenberg model the interchain coupling seems to grow out of super-superexchange via oxygen states.

## Acknowledgments

We are indebted to D. Khomskii and A. Maignan for fruitful discussions. C. Laschinger gratefully acknowledges a Marie Curie fellowship of the European Community program under number HPMT2000-141. This work was supported by the Deutsche Forschungsgemeinschaft (DFG) through Sonderforschungsbereich SFB 484 and by the BMBF (13N6918A).


text:
# Thermal structure of gas in relaxed clusters of galaxies

## 1 Introduction

It has been known for long years that in most clusters of galaxies the intracluster plasma is far from being isothermal, exhibiting strong decrease of the plasma temperature in the vicinity of individual large galaxies. This effect is commonly attributed to ”cooling flows”, that means, increasing radiative energy loss with increasing density of the gas, which is gravitationally attracted by the galaxies. But recent high resolution data from the XMM-Newton and Chandra satellites have shown that this cooling flow model does not fit to the observations. This has newly kindled a vivid discussion on the the origin of the discrepancy.

Already in former years the fact was addressed that thermal conduction would wash out the temperature gradients, if the thermal conductivity of the plasma were of the magnitude given in the classical book by Spitzer (spitzer (1962)). As a possible explanation it has been proposed that magnetic fields may reduce the conduction. Several authors claim very different reduction factors of the thermal conductivity (see e.g.Voigt et al. fab (2002), Markevitch et al. marke (2003)) to fit the data. The existence of magnetic fields in the order of $`0.11\mu \mathrm{G}`$ has been confirmed by various observations (see e.g. Dolag et al. dolag (2001)). But to produce the required reduction, ordered magnetic fields of a still higher magnitude would be necessary, as the heat flux reduction works only perpendicular to the magnetic field lines and there is not yet a generally agreed model to decide, to which extent magnetic fields can influence thermal conductivity. Narayan and Medvedev naray (2001) have shown that fields fluctuating on a large range of length scales can produce only minor changes in the thermal conductivity, while others (see Chandran and Maron chan (2004) and references therein) give reduction factors as high as 10.

One problem with all these reduction factors is, that they would affect the temperature profiles preferably at high temperatures in the outer boundary of the ”cooling flow” region, while radiative cooling is most effective near the centre. But observations exhibit just the opposite effect (Kaastra et al. kaastra (2004))

In all these discussions it is implicitly assumed that in thermodynamic equilibrium conduction leads to a uniform temperature, though it is well known that density gradients may also cause a flux of energy and thus induce temperature gradients, an effect known from the textbooks as ”Dufour effect” (see e.g. Hirschfelder et al. hirsch (1964)). Under terrestrial conditions this effect is regarded mostly as a second order correction to ”normal” conduction, that means, to energy flux due to temperature gradients. We are more familiar with the ”reciprocal” effect, the diffusion of particles caused by temperature gradients, known as thermal diffusion. But under the condition of cosmic plasmas, where stationary density gradients are maintained by gravitational fields, the Dufour effect can considerably alter the equilibrium conditions. Even in relaxed systems there will be stationary temperature gradients to balance the energy flux caused by density gradients.

Normally in a gas we associate ”equilibrium” with constant pressure, density and temperature. But in a system, which is influenced by volume forces such as gravitation, the equilibrium state may well exhibit gradients of the state variables. Mass flow equilibrium in these systems is obtained by the balance between the gravitational force and the counteracting pressure gradient, which is related to the density gradient by some equation of state.

But if there is a density gradient, energy equilibrium can be obtained only, if there is also a temperature gradient, balancing the Dufour effect, to obtain zero net energy flux. Below we will derive the equilibrium conditions for a system, which is fully relaxed with respect to flux of mass, momentum and energy, within the scope of a simplified kinetic model.

In the next section we will discuss the equilibrium conditions in cluster plasmas. In the third section a simplified kinetic model is developed to determine the essential consequences of the Dufour effect. In section 4 the results are compared with observational data from XMM-Newton and Chandra, and in section 5 conclusions are drawn with respect to the ”cooling flow” region as well as for the equilibrium conditions of complete galaxy clusters.

## 2 Equilibrium relations

Speaking of equilibrium means that the net fluxes of all properties determining the state of a system are locally balanced. Complete local equilibrium will of course never be established in a system with radiation losses or with irreversible chemical reactions. But as long as these are not the dominant exchange processes, quasi-equilibrium conditions can be formulated, considering only the dominant transport processes.

To identify the dominant transport processes, it is convenient to compare the time constants for the local change of temperature, density or composition due to some loss and production process or flux. Common to all discussions on the energy balance of the ”cooling flow” region of intracluster gas is, that radiation is regarded as the only energy loss process. This loss has to be balanced by energy gain from inflowing hot matter, by uptake of gravitational energy from inflowing gas, by heat conduction from outer hot regions and energy production by supernovae in galaxies or active galactic nuclei.

The Dufour effect adds a further transport process to the balance, which acts like an additional cooling, opposite to the heating processes mentioned above. It can be regarded as a thermal sedimentation of low energy particles into the vicinity of a central gravitating mass. Microscopically it is closely connected to the ”normal” thermal conduction, as both result from the same process: the transport of kinetic energy of the plasma particles in an inhomogeneous medium. If there is thermal conduction, there is always one fraction caused by temperature gradients and one caused by density gradients. The relation between the two parts will be discussed more in detail in the next section.

To find out, which of the heating or cooling processes preferentially contributes to the balance we compare the time constants of the corresponding processes. In the ”classical” cooling flow model only two processes are considered: radiation losses and inward flow of hot matter. The time constant for radiative cooling is $`\tau _r=(3nkT)/(\mathrm{\Lambda }(T)n^2)`$, where $`\mathrm{\Lambda }(T)`$ is the cooling function (see e.g.Voigt et al. fab (2002)). In the temperature and density range observed in the centre of galaxy clusters typical values are $`kT=3\mathrm{keV}`$ and $`n=10^2\mathrm{cm}^3`$, resulting in a cooling time $`\tau _r=1.2\times 10^{17}\mathrm{sec}`$. As this is much longer than the collapse time $`\tau _c=1/\sqrt{\rho G}=3\times 10^{16}\mathrm{sec}`$, equilibrium can always be reached between gas inflow and radiation loss.

There are, however, two mayor problems with this model. It does not fit to the observations, as the increasing cooling rate should increase with increasing density, so that the temperature profile should exhibit a sharp dip in the centre, which is not observed. Secondly it neglects thermal conduction, which has time constants much shorter than radiation, so that it could wash out all gradients established by radiative cooling. If thermal conduction were at the theoretical value given by Spitzer spitzer (1962) $`\lambda =4.6\times 10^7T^{5/2}\mathrm{erg}\mathrm{cm}^2\mathrm{K}^1\mathrm{sec}^1`$, for a cooling flow region of $`r_c=50\mathrm{kpc}`$ the conduction time scale is $`\tau _\lambda =nkr_c^2/\lambda =10^{16}\mathrm{sec}`$, much less than the cooling time scale. Even if thermal conduction is reduced by the presence of magnetic fields, as proposed by several authors (see Narayan and Medvedev naray (2001), Chandran and Maron chan (2004)), it remains a dominant transport process.

But as has been mentioned above, this does not imply that the temperature should be constant. Due to the opposite effects of temperature and density gradients also under the condition of local thermal equilibrium the temperature will vary in radial direction. In the next section we will derive this equilibrium condition under the assumption that conduction is dominant and radiation losses can be neglected.

## 3 Kinetic model

To determine the equilibrium conditions in a plasma under the influence of an outer force field such as gravitation, we determine the net flux of mass and energy through a control area $`\mathrm{\Delta }F`$ under the condition that the particles in the gas have a Maxwellian distribution everywhere, but with number density and mean energy varying perpendicular to the control area. The net flux is calculated then under the assumption, that any particle crossing the control area retains its velocity and direction, which it has obtained in the last collision, one mean free path from the surface, and in addition is subjected to some acceleration from the force field.

In the energy balance the gain and loss of kinetic energy due to the changing gravitational potential will be omitted, as it cancels out from the balance exactly, when the net particle flux through the control area is zero. As we are interested only in the equilibrium condition, the absolute magnitude of the energy or particle flux is not important, as long as it dominates over other transport, production or loss processes, which appears to be justified for most of the ”cooling flow” region in clusters.

Denoting the direction perpendicular to the surface as z and the angle between this and the direction of particle motion as $`\theta `$, the flux of some property $`G`$ through the surface is

$$j=\mathrm{\Delta }F\mathrm{cos}\theta (j^+j^{})𝑑\mathrm{\Omega },$$
(1)
where $`j^+`$ and $`j^{}`$ are the normal components of the flux of $`G`$ carried by particles moving towards the surface in positive and negative $`z`$ direction under the angle $`\theta `$. The integration has to be taken over half the solid angle $`\pi /2\theta 0`$. Denoting the amount of property $`G`$ carried by particles moving with velocity $`v`$ by $`g(v)`$, the flux components of $`g(v)`$ are

$$j^+=(g(v)\frac{dg(v)}{dz}\lambda \mathrm{cos}\theta )(v\mathrm{cos}\theta +b\frac{\lambda }{v})$$
(2)
$$j^{}=(g(v)+\frac{dg(v)}{dz}\lambda \mathrm{cos}\theta )(v\mathrm{cos}\theta b\frac{\lambda }{v})$$
(3)
$`\lambda `$ is the mean free path, which may depend also on the velocity. The last term denotes the acceleration of the particles during their flight from the last collision to the control surface. Integrating eq.(1) over the solid angle yields the equation

$$j=\pi \mathrm{\Delta }F\left(\lambda v\frac{dg(v)}{dz}\frac{2b\lambda }{v}g(v)\right)$$
(4)
To determine mass and energy transport by particles moving with velocity $`v`$, $`g`$ has to be set to

$$g_\mathrm{M}=n(z)mf(v,z)dv,g_\mathrm{E}=n(z)m\frac{v^2}{2}f(v,z)dv$$
(5)
In case of a neutral gas $`n(z)`$ and $`m`$ are the numerical density and mass of atoms, $`f(v,z)`$ is the distribution function

$$f(v)dv=\frac{4}{\sqrt{\pi }}\frac{v^2}{\beta ^3}e^{\frac{v^2}{\beta ^2}}dv$$
(6)
$`\beta (z)`$ being the abbreviation $`\beta (z)=\sqrt{2kT(z)/m}`$. To calculate the quantity $`dg/dz`$ we also need the derivative

$$\frac{df}{dz}=\frac{d\beta }{dz}\frac{4}{\sqrt{\pi }}\left(\frac{2v^4}{\beta ^6}\frac{3v^2}{\beta ^4}\right)e^{\frac{v^2}{\beta ^2}}=\frac{1}{\beta }\frac{d\beta }{dz}\left(\frac{2v^2}{\beta ^2}3\right)f(v)$$
(7)
While in a neutral gas the parameter $`b`$ in eq.(4) is the normal gravitational acceleration, in a plasma the action of the force is somewhat indirect. Energy transport is mediated preferably by electrons. Thus in eqs.(5) $`n(z)`$ and $`m`$ are the electron density and mass. But in this case acceleration by the external force is acting also on the ions, which transfer the force to the electrons by Coulomb interaction. Thus the parameter $`b`$ in eqs.(2) and (3) is not the gravitational acceleration of free moving electrons. But for the equilibrium condition of zero transport the absolute magnitude of $`b`$ is not relevant, when we only want to know, under which conditions the net flux of electrons and the flux of energy carried by them are in balance.

An additional difference between neutral gas and plasma results from the fact that the collision cross section between neutral atoms is nearly independent of the collision energy, so that the mean free path $`\lambda `$ does not depend on $`v`$. In plasmas, where Coulomb interaction is dominant, the mean free path increases with $`v^2`$. Thus in the further calculations we set $`\lambda =\lambda _0(v/v_0)^\alpha `$ with $`\alpha =0`$ for neutral gas and $`\alpha =2`$ for a fully ionised plasma.

Introducing eqs.(5) to (7) into eq.(4), we get for the mass and energy flux

$$j_\mathrm{M}=2\beta ^\alpha \left(x\beta \frac{dn}{dz}+x(2x^23)n\frac{d\beta }{dz}\frac{2bn}{\beta }\frac{1}{x}\right)h(x)dx$$
(8)
$$j_\mathrm{E}=\beta ^{\alpha +2}\left(x^3\beta \frac{dn}{dz}+x^3(2x^23)n\frac{d\beta }{dz}\frac{2bn}{\beta }x\right)h(x)dx$$
(9)
with $`x=v/\beta `$ and $`h(x)=2\sqrt{\pi }m\mathrm{\Delta }F(\lambda _0/v_0^\alpha )x^{2+\alpha }e^{x^2}`$.
Integrating these equations over all x, we obtain the conditions for zero mass and energy flux:

$$\beta \frac{dn}{dz}F_{3+\alpha }+n\frac{d\beta }{dz}[2F_{5+\alpha }3F_{3+\alpha }]\frac{2bn}{\beta }F_{1+\alpha }=0$$
(10)
$$\beta \frac{dn}{dz}F_{5+\alpha }+n\frac{d\beta }{dz}[2F_{7+\alpha }3F_{5+\alpha }]\frac{2bn}{\beta }F_{3+\alpha }=0$$
(11)
The abbreviation $`F_k`$ stands for the integral $`F_k=_0^{\mathrm{}}x^ke^{x^2}𝑑x`$. The different arguments of $`F_k`$ in eqs.(10) and (11) can be eliminated by the recurrence formula $`F_{k+2}=(k+1)/2\times F_k`$, so that we finally obtain the equilibrium conditions for mass and energy transport:

$$\frac{2+\alpha }{2}\beta \frac{dn}{dz}+\frac{(2+\alpha )(1+\alpha )}{2}n\frac{d\beta }{dz}\frac{2bn}{\beta }=0$$
(12)
$$\frac{4+\alpha }{2}\beta \frac{dn}{dz}+\frac{(4+\alpha )(3+\alpha )}{2}n\frac{d\beta }{dz}\frac{2bn}{\beta }=0$$
(13)
It is immediately evident that this is a set of linear equations for the gradients $`d\beta /dz`$ and $`dn/dz`$. If there are no external forces $`(b=0)`$, there exists only the trivial solution $`d\beta /dz=0`$ and $`dn/dz=0`$.

In the presence of volume forces such as gravitation eqs.(12) and (13) are inhomogeneous and allow a non-trivial solution, which constitutes a fixed relation between temperature gradient and particle density gradient or, because of quasi-neutrality, also between temperature and mass density gradient. Eliminating $`b`$ leads to the relation

$$\beta \frac{dn}{dz}+(2\alpha +5)n\frac{d\beta }{dz}=0$$
(14)
with the solution $`n\beta ^{(2\alpha +5)}=\mathrm{const}.`$ Replacing the number density by the mass density $`\rho `$ and $`\beta `$ by $`\sqrt{2kT/m}`$ we finally find the condition $`\rho T^{5/2}=\mathrm{const}.`$ for neutral gas and $`\rho T^{9/2}=\mathrm{const}.`$ for a fully ionised plasma.

This result is independent from the actual magnitude of the coefficient of thermal conduction, if it is influenced by magnetic fields or not. The only conditions, which have to be satisfied, are the existence of local thermal equilibrium and the condition, that transport of kinetic energy of particles is the dominant transfer mechanism.

## 4 Observations

During the last years several spatially resolved x-ray emission measurements of galaxy clusters have become available, especially from the Chandra and XMM-Newton satellites. Though temperature determinations are still rather poor due to projection effects and necessary background corrections, intensity profiles and from these density profiles have been obtained with high accuracy. But in spite of the remaining error margin the classical cooling flow model can definitely be ruled out as an explanation of the observed data (see Kaastra et al. kaastra (2004)). Neither the radial temperature profiles fit to those expected from this model nor do the spectra show the expected intensity of metal lines of lower ionisation levels.

Good agreement with observations can, however be obtained with a model based on thermal conduction with inclusion of the Dufour effect. We have applied the model to the data given by Kaastra et al. kaastra (2004). In fig.1 the quantity $`T\times \rho ^{2/9}`$ is plotted for the seven most luminous clusters, for which data are given by Kaastra et al. Only the clusters with the highest X-ray flux ($`>5\times 10^{11}\mathrm{erg}\mathrm{cm}^2\mathrm{sec}^1`$, 0.1-2.5 keV) have been used, as for these clusters the temperature profiles appear sufficiently accurate. The data points in fig.1 refer to temperatures in consecutive annular channels starting from centre to the boundary of the cooling region. (For details see Kaastra kaastra (2004)). For all seven clusters the constancy of $`T\times \rho ^{2/9}`$ over the cooling region is evident. This can be regarded as a confirmation that thermal conduction with inclusion of the Dufour effect is the dominant energy transfer mechanism and that radiative cooling is less important.

Data from Perseus cluster have also been obtained by Chandra (see Sanders et al. sanders (2004)). Fig.2 shows a comparison of the temperature profiles obtained from direct measurements by Chandra and XMM-Newton compared with data derived from the corresponding density profiles, applying $`\rho T^{9/2}=\mathrm{const}.`$ The curves are fitted at the cluster centre. As can be seen, the curves derived from density data agree very well with each other and with the direct measurements from XMM-Newton, while the temperatures from Chandra deviate considerably in the outer region. It should be mentioned, however, that the Chandra data are obtained by averaging angular resolved profiles which exhibit a large spread ($`>1.5`$ keV) in the outer region.

Whether radiative loss can completely neglected or has at least to be added as a correction to the temperature profile resulting from pure heat conduction, cannot be decided from the present data. For the model this would require a better knowledge of the absolute magnitude of the thermal conductivity. This requires an accurate theory of the reduction effects by magnetic fields and more knowledge on the distribution of these fields over the cluster volume.

## 5 Conclusions

We have shown that thermal conduction appears to be the dominant energy transport mechanism and that a pure conduction model can describe the observed temperature profiles, when energy flux induced by density gradients is added to the balance. The relative temperature gradients are small compared to the associated density gradients. But temperature changes by a factor of two or three, as they have been found in the central regions of the clusters mentioned above and which are reported also for other ”cooling flow” clusters by several authors (see e.g. Voigt et al. fab (2002)), can well be understood from the Dufour effect, as the gas density may change by a few orders of magnitude in the vicinity of large galaxies. It appears likely that most of the observed temperature drop in the ”cooling flows” around these galaxies is not caused by radiative cooling but by the ”conductive cooling” associated with the density gradients.

It should be noted that the Dufour effect is not restricted to the cores of galaxies, but it should be important also in the hot corona of individual stars, where it can help to explain the observed temperature gradient from the surface into the interstellar medium.

Applying the model to complete galaxy clusters leads to the result that the plasma temperature should continuously rise with increasing distance from the centre. Observational data are contradictory in the outer region. Some authors claim a temperature drop at the outskirts of the observable region, others find constant or increasing temperature (see Mushotzky mush (2004) and references given there). Most of the recent papers find increasing temperatures at least out to the point, where the increase of the error bars exceedes the expected temperature gradients.

Also the conduction model should be taken with care at the low densities in the outer regions, as the conditions of local thermal equilibrium are no longer satisfied due to the extremely large mean free path of electrons. But it can be stated that, if thermal conductivity were close to the Spitzer value, clusters would have cooled off in the time since formation, unless the clusters are embedded in a very hot but thin intergalactic plasma. The existence of a diffuse high energy x-ray background would be a strong hint to the existence of such an intergalactic plasma.

Measurements of the background in the energy range above that of the intracluster gas are rare. Observations from the HEAO-1 satellite (Boldt boldt (1987)) have found a continuous background, which can be interpreted as bremsstrahlung in the temperature range of 40 to 50 keV, but the spatial resolution of these early data is so poor that it cannot be decided, if it is really radiation from a continuously distributed plasma or from the superposition of individual point sources. As a plausible explanation for the formation of a hot continuous background plasma one could think of the relativistic matter jets expelled by quasars. During the time period of several Gyr this matter should have become homogeneously distributed in space.

It has been shown that the conduction model presented in this paper can explain some observations in the central regions of the intracluster plasma, but it deserves more confirmation by detailed measurements, until the astrophysics community will be convinced. Also the extension to the outer regions of galaxy clusters would bring additional confirmation, but this requires the extension of spatially resolved measurements into the high energy range, as they are planned with the launch of future satellites like Astro-E2, which are scheduled for the next years.