text:
# Electronic Structure of Calcium Hexaboride within the Weighted Density Approximation

## Abstract

We report calculations of the electronic structure of CaB6 using the weighted density approximation (WDA) to density functional theory. We find a semiconducting band structure with a sizable gap, in contrast to local density approximation (LDA) results, but in accord with recent experimental data. In particular, we find an $`X`$-point band gap of 0.8 eV. The WDA correction of the LDA error in describing the electronic structure of CaB6 is discussed in terms of the orbital character of the bands and the better cancelation of self-interactions within the WDA.

The discovery of ferromagnetism in the lightly doped CaB6 and SrB6 resulted in considerable experimental and theoretical work aimed at understanding this unexpected phenomenon. The reported ferromagnetism is particularly remarkable considering that the material, Ca1-xLaxB6, contains no elements with partially filled $`d`$ or $`f`$ orbitals, and even more strangely, the ferromagnetism is associated with very low doping levels ($`<1\%`$), and therefore low carrier concentrations, has a very low moment ($``$ 0.07 $`\mu _\mathrm{B}`$/La for $`x`$=0.005), but a very high Curie temperature, $`T_C`$ as high as 600K (Ref. young99, ). Further, a semiconducting material, with ferromagnetism of this type, could potentially be of importance in producing room temperature spintronic devices. wolf ; tromp2001

Two different classes of theories were advanced early on in an effort to explain the strange ferromagnetism. The first was based on the doping of an excitonic or like system, i.e. an unconventional state in the parent, CaB6, zhit99 ; bal2000 ; bar2000 ; mura2002  and the second was based on the ferromagnetic phase of the dilute electron gas. young99 ; cep99 ; ortiz99  In the first class, the ferromagnetism arises from La induced doping of carriers into an excitonic ground state. This depends on a bare band structure that is semimetallic or very close to it, so that the prerequisite instability of CaB6 to an excitonic state can exist. The second class mentioned does not depend on such an electronic structure, but rather involves a ferromagnetic instability of a low density electron gas, though this class of theories suffers from difficulties obtaining the very high reported values of $`T_C`$. More pedestrian explanations for the observed ferromagnetism involving defects and their interactions and also extrinsic effects have been advanced more recently. gavilano ; matsubayashi ; matsubayashi2 ; young2

This naturally resulted in a strong interest in the electronic structure of the parent compound CaB6. First-principles calculations within the local density approximation (LDA) predicted that CaB6 and SrB6 are semimetals with a small overlap at X point of the simple cubic Brillouin zone, or very small gap insulators, mass97 ; rod2000  consistent with the excitonic theories being formulated at that time. However, it is well known that the LDA often underestimates band gaps, and in some extreme cases it can predict a semimetallic band structure for a semiconductor, one example being Ge. fil  In simple semiconductors, the LDA underestimates of the band gaps are understood as being associated with an extreme non-analytic and non-local behavior of density functional theory as the particle number is changed to include a single conduction electron, godby ; hybertsen  implying the need for a self-energy correction, or at least an orbital dependent potential to obtain a realistic description of band gaps.

The many-body quasiparticle $`GW`$ approximation hedin  is very succesful in describing the band structures of insulators. In particular the self-energy correction implicit in it corrects most of the band gap errors in simple semiconductors and other materials as well. aryasetiawan  However, its results for CaB6 have differed depending on the formulation of $`GW`$ being used. Pseudopotential based $`GW`$ calculations tromp2001  yielded a sizable 0.8 eV band gap at X point, while recent all-electron $`GW`$ results kino2002  showed an unexpected increase in band overlap relative to the LDA, perhaps related to the treatment of deep-core states. kotani2001 ; ku2002

The electronic structure of CaB6 has also been the subject of some debate from an experimental point of view, although a semiconducting electronic structure has become more generally accepted now. de Haas-van Alphen (dHvA) good98 ; hall2001  and Shubnikov de Hass (SdH) aron99  experiments were consistent with a semimetallic nature, while recent angle-resolved photoemission spectroscopy (ARPES) and resonant inelastic x-ray scattering (RIXS) measurements imply that CaB6 is a semiconductor with a 1.0 (Ref. souma2003, ) - 1.15 (Ref. denlinger2002, ) eV band gap. Recently, Rhyee et al. reported transport, optical, and tunneling measurements all consistent with a 1.0 eV band gap for pure CaB6. They argued that the previously measured semimetallic characteristics in CaB6, SrB6, and EuB6 actually originate from boron vacancies. rhyee2003

We note that, as pointed out by Rodriguez and co-workers, rod2000 ; pic03  the valence and conduction bands of CaB6 and the related hexaborides are of very different character, namely B $`p`$ valence bands and Ca derived conduction bands. This is reminiscent of the situation in the rare-earth trihydrides, e.g. YH3. Those materials feature a semiconducting ground state, which is however predicted to be metallic in the LDA due to a band overlap between nominally occupied H derived states and nominally unoccupied metal states. In YH3, this error is understood to be due to incomplete cancelation of self-interaction errors on the H sites in the LDA. These self-interactions lead to a H band position that is too high relative to the conduction bands, which is at least partly corrected in $`GW`$ calculations. eder ; vg2 ; miyake  It should be noted though that such self-interactions are a local effect, more akin in some ways to the on-site Coulomb effects that lead to Mott insulating states in transition metal oxides, than to the non-local, non-analytic discontinuities in density functional theory that are important in simple semiconductors. cohen  As such, the rare earth trihydrides and perhaps CaB6 are interesting test cases for theories of correlated electron systems, including density functional based approaches.

Here, we present results of first-principles calculations based on the density functional theory (DFT) within a non-local density approximation, namely the weighted density approximation (WDA). gunna76 ; alonso78 ; gunna79 ; singh93  The WDA uses a model pair-distribution function to incorporate non-local information about the charge density into the exchange correlation energy and potential, via an exact expression for the total energy. Unlike $`GW`$ and like approaches, the WDA is oriented towards total energies rather than excitation spectra. As such, not surprisingly, it improves ground state properties, like lattice parameters, relative to the LDA. Like the LDA, the WDA is exact for the total energy of the uniform electron gas, but unlike the LDA it is also exact for single electron systems, where it is self-interaction free. We emphasize that the WDA is a smooth non-pathalogical approximation to density functional theory, with only short range non-locality. The WDA exchange correlation potential is a standard orbital independent (for valence states) potential, which different orbitals feel differently only through differences in the spatial distrubutions of their charge densities. It does not include the discontinuities that are present in the exact density functional theory and are responsible for most of the band gap correction in simple semiconductors like Si. Detailed calculations for YH3 and LaH3 have shown that the WDA can yield insulating gaps of similar magnitude to the $`GW`$ approximation in these hydrides due to the better cancelation of self-interaction errors, while at the same time improving the quality of ground state properties. wu2004

In the WDA the exchange-correlation (xc) energy expression is very similar to the general DFT form,

$$E_{\mathrm{xc}}^{\mathrm{WDA}}[n]=\frac{n(𝐫)n(𝐫^{})}{|𝐫𝐫^{}|}G[|𝐫𝐫^{}|,\overline{n}(𝐫)]𝑑𝐫𝑑𝐫^{},$$
(1)
where $`G`$ is the approximate model coupling constant averaged pair-distribution function, and $`\overline{n}(𝐫)`$ is the weighted density which can be determined from the sum rule which states that the charge density in the exchange correlation hole is unity.

$$n(𝐫^{})G[|𝐫𝐫^{}|,\overline{n}(𝐫)]𝑑𝐫^{}=1.$$
(2)
The WDA is exact for both a one-electron system and the uniform electron gas, independent of the model $`G`$ used, provided that it is monotonic and properly normalized. Because the exact $`G`$ function is not known, different forms of $`G`$ have been tried, mazin  from which we chose the following one

$$G(r,n)=c\{1\mathrm{exp}([\frac{r}{\lambda }]^4)\},$$
(3)
This choice yields good energetics (structural properties) and band structures for the rare earth trihydrides. wu2004  This type of $`G`$ was derived from the Gunnarsson-Jones gunna80  ansatz, but it has a shorter range than the original one. Details of the WDA formalism used in this paper can be found in Refs. wu2004,  and wu2003, .

The calculations were done using a planewave basis with pseudopotentials. The k-grids for the Brillouin zone samplings during iteration to self-consistency were as dense as $`12\times 12\times 12`$, and the energy cutoff was 121 Ry. dodpw  The Ca semicore states were treated as valence states for the purpose of constructing Troullier-Martins pseudopotentials. Shell partitioning was employed as discussed in Ref. singh93, . The crystallographic structure of CaB6 is simple cubic, with Ca atoms located at cubic corners. The B atoms form octahedral cages around the cube center. The internal parameter $`x`$ determines the distance between a B atom and the cubic center. When $`x=0.207`$, the intra- and inter- octahedral B-B distances are equal. Our calculated lattice constant $`a`$ and internal parameter $`x`$ within the LDA and WDA are shown in Table 1. The WDA results are in excellent agreement with experimental data, while the LDA predicts a lattice constant 1.6% smaller, which is a typical LDA tendency. Thus, reflecting its origins in density functional theory, where energy and density are primary variables, the WDA improves upon the LDA structure.

To facilitate comparison of the LDA and WDA electronic band structures, we used the experimental values of $`a`$ and $`x`$ for both the LDA and WDA. Fig. 1 shows our results for the band structures of CaB6. The LDA (left panel) predicts a tiny band gap (slightly more than 0.1 eV) at the $`X`$-point, which agrees with Ref. kino2002 , while other LDA results show a tiny overlap mass97 ; rod2000 ; tromp2001 . This small difference in band structure is due to different choices of the LDA scheme, chemistry (Sr vs. Ca) and lattice structure as well as the approximations in using pseudopotentials. To check this we did well converged general potential linearized augmented planewave (LAPW) calculations within the LDA for two structures. The first is the same experimental structure as used for Fig. 1, and the second is a structure where we relaxed the internal coordinate $`x`$ in the LDA keeping the lattice parameter fixed at the experimental value. This yielded $`x`$=0.2025. With the experimental value of $`x`$=0.207 the LAPW calculations yielded a small band gap of 0.2 eV, while at $`x`$=0.2025 a semimetal with a 0.1 eV band overlap is predicted. This implies that the structure is the key ingredient in reconciling the various results and also gives an error estimate of $``$ 0.1 eV for our pseudopotential calculations. On the other hand, the right panel shows that within the WDA, the bands are shifted, quite similar to the cases of YH3 and LaH3, resulting in a sizable band gap of 0.8 eV at X point. Parallel WDA calculations using the uniform electron gas $`G`$ function also yield a gap opening, though by a smaller amount. The band gap in that case is 0.3 eV.

Thus the WDA, with an appropriate choice of $`G`$, predicts CaB6 to be a semiconductor with a gap comparable to that found in standard $`GW`$ calculations tromp2001  and in experimental spectroscopic measurements. souma2003 ; denlinger2002 ; rhyee2003  This shows that there exists a physical approximation to density functional theory that is continuous and has only short range non-locality and that gives both a reasonable band gap and ground state properties (crystal structure) in agreement with experiment for CaB6. We expect that WDA calculations would yield band gap opening in SrB6 as well.

In summary, WDA calculations predict that CaB6 is a semiconductor, and the calculated band gap is in excellent agreement with recent experiments and quasi-particle $`GW`$ calculations. Thus the WDA improves both the structural properties and electronic structure of these hexaborides, implying that the main error in the LDA description is associated with local self-interaction effects on B rather than the neglect of the non-analytic discontinuity in density functional theory across the band gap.

We are grateful for discussions with M.H. Cohen, M. Gupta, R. Gupta, I.I. Mazin and W.E. Pickett. DJS is grateful for the hospitality of the University of Paris-Sud where some of the ideas leading to this work were developed. Work at the Geophysical Laboratory was supported by Office of Naval Research (ONR) grants N000149710052, N00014-02-1-0506 and N0001403WX20028. Work at the Naval Research Laboratory is supported by ONR. Calculations were done on the Center for Piezoelectrics by Design (CPD) computer facilities.


format: text

text:
# Stellar and gaseous abundances in M82Based on observations with the Near Infrared Camera Spectrometer (NICS) mounted at the Telescopio Nazionale Galileo (TNG), La Palma, Spain and archival XMM-Newton and Chandra data.

## 1 Introduction

The signature of the star formation (SF) history of a galaxy is imprinted in the abundance patterns of its stars and gas. Determining the abundance of key elements released in the interstellar medium (ISM) by stars with different mass progenitors and hence on different time scales, will thus have a strong astrophysical impact in drawing the global picture of galaxy formation and evolution (McWilliam, 1997).

So far, the metallicities of distant starburst (SB) galaxies have been mainly derived by measuring the nebular lines associated to their giant HII regions (Storchi-Bergmann, Calzetti & Kinney, 1994; Coziol et al., 1999; Henry, Edmunds & Köppen, 2000). Direct determination of O and N abundances can be obtained once electron temperature (Te) is well established (McCall, 1984). At high metallicities Te decreases due to strong lines cooling and its major diagnostic (the weak \[OIII\] $`\lambda `$4636) is no longer observable. Several alternative methods based on other strong lines were proposed in the past years (Pagel et al., 1979; Dopita & Evans, 1986; McGaugh, 1991) to derive Te, but all of them rely on a large number of physical parameters and model assumptions for the geometry and the nature of the ionizing populations (c.f. e.g Stasinska (2001) for a review). Depending on the adopted calibrations, quite different abundance sets can reproduce the observed line ratios.

However, abundances in SB galaxies can be also measured from absorption stellar features in the near IR and/or emission lines in the hot X-ray gas. Stellar abundances can be poorly constrained from optical spectra, since the nebular emission strongly dilutes the absorption lines and dust can heavily obscure the central regions where most of the burst activity is concentrated. On the contrary, in the near IR the stellar continuum due to red supergiants (RSGs) usually largely dominates over the possible gas and dust emission (Oliva & Origlia, 1998; Origlia & Oliva, 2000). The X-ray emission from SB galaxies at E$`>`$2 keV is mainly due to high mass X-ray binaries with power-law spectra, while at lower energies hot plasma mainly heated by SN explosions (Dahlem et al., 1998) dominates and the resulting spectrum is a superposition of thermal continuum and emission lines from atomic recombinations (Fabbiano, 1989; Persic & Rephaeli, 2002).

Determining the abundances of SB galaxies is not only relevant to constrain their SF history, but also offers the unique chance of directly witnessing the enrichment of the ISM (see e.g. Maeder & Conti, 1994, for a review). Metals locked into stars give a picture of the enrichment just prior to the last burst of SF, while the hot gas heated by SNe II explosions should trace the enrichment by the new generation of stars. Nebular metallicities from the cold gas should potentially trace the stellar ones if the cooling and mixing timescales are slow (up to 1 Gyr), as suggested by some models (see e.g. Tenorio-Tagle, 1996) since the gas enriched by the ongoing burst is still too hot to be detected in the optical lines. On the contrary, if rapid (a few Myr) cooling and mixing occur as suggested by other models (see e.g. Recchi, Matteucci & D’Ercole, 2001), the nebular abundances could trace the gas enriched by the new population of stars (similarly to X-rays). Hence, important constraints on the cooling and mixing timescales of the gas can be provided by comparing the metallicity inferred from stars and hot plasma with those of the cold gas.

Sect. 2 briefly reviews the nebular abundance estimates of M82, as published so far in the literature. Sect. 3 shows our near IR spectra of the nuclear region of M82 and the derived stellar abundances. Sect. 4 presents a re-analysis of the XMM-Newton and Chandra nuclear X-ray spectra and the hot gas–phase abundances. Sect. 5 summarized the overall chemical enrichment scenario in M82 as traced by the stellar, cold and hot gas-phase abundances, while in Sect. 6 we draw our conclusions.

## 2 Nebular abundances of M82 from the literature

M82 is a prototype of SB galaxies (Rieke et al., 1980), experiencing a major SF formation episode in its nuclear regions, with strong super-wind and SN activity (Wills et al., 1999). Despite the many multi-wavelength studies (see e.g. Golla, Allen & Kronberg, 1996, and references therein), only few measurements of nebular abundances, mainly derived from mid/far IR forbidden lines, have been published so far. From the analysis of the \[SiII\] 35$`\mu `$m line in the central 300 pc in radius at the distance of M82 (i.e. d$``$3.6 Mpc), Lord et al. (2001) suggest an indicative $`[\mathrm{Si}/\mathrm{H}]0.4`$ dex relative to the Solar value of Grevesse & Sauval (1998).

More recently, Förster Schreiber et al. (2001), by analyzing the nuclear spectra of M82 taken with the Short Wavelength Spectrometer (SWS) onboard the Infrared Space Observatory (ISO) within an aperture corresponding to $`100\times 200`$ pc on the sky, give cold gas-phase abundances of \[Ne/H\]=+0.1, \[Ar/H\]=+0.3 and \[S/H\]=–0.7 dex relative to the Solar values of Grevesse & Sauval (1998), with an overall uncertainty between 0.2 and 0.3 dex.

By taking the line ratios listed by Alloin, Collin-Souffrin, & Joly (1979) which refer to the central $``$100 pc in radius of M82 (Peimbert & Spinrad, 1970), and the empirical calibrations of McGaugh (1991) which give N and Fe abundances as a function of \[O/H\] (see e.g. their Figure 11, equations 9 and 10a), one can also obtain rough estimates of \[O/H\]=+0.1, \[N/H\]=0.0 and \[Fe/H\]=–0.3 dex relative to the Solar values of Grevesse & Sauval (1998), with an overall uncertainty of $`\pm `$0.3 dex.

The first attempt to measure hot gas-phase abundances in M82 has been made quite recently by Ptak et al. (1997) and slightly revised by Umeda et al. (2002), by analyzing ASCA X-ray (0.4-10.0 keV) spectra. Fe and O abundances as low as 1/20 Solar and significantly higher (by a factor between 3 and 10) Si, S, Mg, Ca ones have been obtained, also consistent with BeppoSAX results (Cappi et al., 1999).

More recently, Read & Stevens (2002) by analyzing the spectra obtained with the Reflection Grating Spectrometer (RGS) onboard XMM-Newton found near-Solar Fe and O and super-Solar (between a factor of 2 and 5) Mg, Si, Ne, N abundances. While the abundance ratios are in reasonable agreement within $`\pm `$0.3 dex with Ptak et al. (1997), there is one order magnitude difference between the zero-point of the two calibrations.

## 3 IR spectra and stellar abundances

Near-IR spectroscopy is a fundamental tool to obtain accurate abundances of key elements like Fe, C, O and other metals (e.g. Si, Mg, Ca, Na, Al, Ti) in cool stars ($`\mathrm{T}_{\mathrm{eff}}5000`$ K). Several atomic and molecular lines are strong, not affected by severe blending, hence also measurable at low-medium resolution (Kleinmann & Hall, 1986; Origlia et al., 1993; Wallace & Hinkle, 1997; Joyce et al., 1998; Frogel et al., 2001), making them powerful abundance tracers not only in stars but also in more distant stellar clusters and galaxies and in a wide range of metallicities and ages (Origlia et al., 1997; Oliva & Origlia, 1998). Abundance analysis from integrated spectra of galaxies requires full spectral synthesis techniques to properly account for line blending and population synthesis to define the dominant contribution to the stellar luminosity.

In SB galaxies the stellar IR continuum usually dominates over the possible gas and dust emission. This represents a major, conceptual simplification in population and spectral synthesis techniques, making possible and easier the interpretation of integrated spectra from distant stellar clusters and galaxies.

Since the near IR continuum of SB galaxies is dominated by luminous RSGs (see e.g. Origlia & Oliva, 2000, for a review), their integrated spectrum can be modeled with an equivalent, average star, whose stellar parameters (temperature Teff, gravity logg and microturbulence velocity $`\xi `$) mainly depend on the stellar age and metallicity. Both observations and evolutionary models (see e.g. Keller, 1999; Origlia et al., 1999; Origlia, 2003, and reference therein) suggest that RSGs of ages between $``$6 and 100 Myr and metallicities between 1/10 and Solar are characterized by low gravities (log g$`<`$1.0), low temperatures ($``$ 4000 K) and relatively high microturbulence velocity ($`\xi `$3 km/s). A variation in the adopted stellar parameters for the average RSG population by $`\mathrm{\Delta }`$Teff$`\pm `$200 K, $`\mathrm{\Delta }`$log g=$`\pm `$0.5 and $`\mathrm{\Delta }\xi `$$``$0.5 km s-1 implies a $`\pm 0.1`$ dex change in the abundances estimated from atomic lines, and $`\pm 0.2`$ dex in those estimated from the molecular lines, which are more sensitive to stellar parameters.

1.0–1.8$`\mu `$ longslit spectra (see Fig. 1) at R=2,500 of the nuclear region of M82 were obtained with the Near IR Camera Spectrometer (NICS) mounted at the Italian Telescopio Nazionale Galileo (TNG) on December 2002. The spectra were sky-subtracted, flat-fielded and corrected for atmospheric absorption using an O-star spectrum as reference. They were wavelength calibrated by using a Ne-Ar lamp and the monodimensional spectra were extracted by summing over the central $`0.5\mathrm{}\times 3.0\mathrm{}`$, corresponding to an aperture projected on the source of $`9\times 52`$ pc at the distance of M82. The spectra have been normalized to the continuum, which was determined applying a low-pass smoothing filter to each spectrum. Total integration times of 72 and 32 min in the J and H bands, respectively, were used, providing a final signal to noise ratio $``$40.

By measuring the absorption line broadening a stellar velocity dispersion $`\sigma 105\pm `$20 km/s has been derived, in perfect agreement with the estimate by Gafney, Lester & Telesco (1993) and in the typical range of values measured in other massive SB galaxies (Oliva et al., 1999).

A grid of synthetic spectra of red supergiant stars for different input atmospheric parameters and abundances have been computed, using an updated (Origlia et al., 2002, 2003) version of the code described in Origlia et al. (1993). Briefly, the code uses the LTE approximation and is based on the molecular blanketed model atmospheres of Johnson, Bernat & Krupp (1980). It includes several thousands of near IR atomic lines and molecular roto-vibrational transitions due to CO, OH and CN. Three main compilations of atomic oscillator strengths are used, namely the Kurucz’s database (c.f. http://cfa-www.harward.edu/amdata/ampdata/kurucz23/sekur.html), and those published by Bièmont & Grevesse (1973) and Meléndez & Barbuy (1999).

Our code provides full spectral synthesis over the 1-2.5 $`\mu `$m range. Given the high degree of line blending, our abundance estimates are mainly obtained by best-fitting the full observed spectrum and by measuring the equivalent widths of a few selected features (cf. Fig. 1), dominated by a specific chemical element, as a further cross-check. The equivalent widths have been measured by performing a Gaussian fit with $`\sigma `$ equal to the measured stellar velocity dispersion, typical values ranging between 0.5 and 3 Å, with a conservative error of $`\pm `$200 mÅ to also account for a $`\pm `$2% uncertainty in the continuum positioning.

By fitting the full observed IR spectrum and by measuring the equivalent widths of selected lines, we obtained the following best-fit stellar parameters and abundance patterns for M82: Teff=4000, logg=0.5, $`\xi `$=3, \[Fe/H\]=–0.34; \[O/Fe\]=+0.34, \[$`<`$Si,Mg,Ca$`>`$/Fe\]=+0.38; \[Al/Fe\]=+0.5; \[C/Fe\]=–0.26; 12C/13C$`<`$10. Table 1 lists the derived abundances and their associated random errors at 90% confidence. Reference Solar abundances are from Grevesse & Sauval (1998).

Synthetic spectra with lower element abundances are systematically shallower than the best-fit solution, while the opposite occurs when higher abundances are adopted. In order to check the statistical significance of our best-fit solution, as a function of merit we adopt the difference between the model and the observed spectrum (hereafter $`\delta `$). In order to quantify systematic discrepancies, this parameter is more powerful than the classical $`\chi ^2`$ test, which is instead equally sensitive to random and systematic scatters (see Origlia et al., 2003, for more details).

Since $`\delta `$ is expected to follow a Gaussian distribution, we compute $`\overline{\delta }`$ and the corresponding standard deviation ($`\sigma `$) for the best-fit solution and 6 test models with abundance variations $`\mathrm{\Delta }[X/H]=\pm `$0.1, 0.2 and 0.3 dex with respect to the best-fit. We then extract 10000 random subsamples from each test model (assuming a Gaussian distribution) and we compute the probability $`P`$ that a random realization of the data-points around a test model display a $`\overline{\delta }`$ that is compatible with the best-fit model. $`P1`$ indicates that the model is a good representation of the observed spectrum.

The left panel of Fig. 2 shows the average results for the observed J and H band spectra of M82. It can be easily appreciated that the best-fit solution provides in all cases a clear maximum in $`P`$ ($`>`$99%) with respect to the test models. More relevant, test models with an abundance variation $`\mathrm{\Delta }[X/H]\pm 0.2`$ dex lie at $`1.5\sigma `$ from the best-fit solution, while test models with $`\mathrm{\Delta }[X/H]\pm 0.3`$ dex lie at $`\sigma 2`$ from the best-fit solution.

The analysis of the line equivalent widths provide fully consistent results. The right panel of Fig. 2 shows the average difference between the model and the observed equivalent width measurements. Models with $`\pm 0.2`$ dex abundance variations from the best-fit solution are still acceptable at a $`1.5\sigma `$ significance level, while those with $`\pm 0.3`$ dex variations are only marginally acceptable at a $`23\sigma `$ level.

Models with stellar parameters varying by $`\mathrm{\Delta }`$Teff$`\pm `$200 K, $`\mathrm{\Delta }`$log g=$`\pm `$0.5 and $`\mathrm{\Delta }\xi `$$``$0.5 km s-1 and abundances varying accordingly by 0.1-0.2 dex, in order to still reproduce the deepness of the observed features, are also less statistical significant (on average only at $`2\sigma `$ level) with respect to the best-fit solution. Hence, as a conservative estimate of the systematic error in the derived best-fit abundances, due to the residual uncertainty in the adopted stellar parameters, one can assume a value of $`\pm 0.1`$ dex.

By taking into account the overall uncertainty in the definition of the average population and the statistical significance of our spectral synthesis procedure, we can safely conclude that the stellar abundances can be constrained well within $`\pm `$0.3 dex and their abundance ratios down to $``$0.2 dex, since some (if not all) of the stellar parameter degeneracy is removed.

## 4 X-ray spectra and hot gas-phase abundances

The determination of the abundances of the hot X-ray emitting gas in SB galaxies has traditionally suffered from large uncertainties. Indeed, the low angular and spectral resolution of the various X-ray telescopes in the pre-XMM-Newton/Chandra era did not allow to disentangle between point sources and hot gas emission, making abundance determinations severely model-dependent (Dahlem et al., 2000). The problem lies in the fact that the X–ray spectrum of SB galaxies contains two major components: the emission from hot diffuse gas and the integrated contribution of point sources. The first is described by an optically thin thermal spectrum plus emission lines, the latter by a power-law. If the angular resolution is not good enough, it is not possible to reliably subtract the point sources from the total spectrum, so that the equivalent widths of the emission lines (and thus the element abundances) are not unambiguously defined.

On the other hand, high angular resolution alone is not enough: recent Chandra  studies (Strickland et al. 2002; Martin, Kobulnicky & Heckman 2002) demonstrated that the relatively low spectral resolution of the ACIS detector makes individual element abundance analysis still problematic. The high spectral resolution of the Chandra gratings rapidly degrades for sources which are much more extended than the instrumental Point Spread Function (PSF) and thus these gratings are almost useless for the study of nearby SB galaxies. However, by coupling the high angular resolution ($`0.5^{\mathrm{}}`$ PSF) of Chandra with the high spectral resolution ($`\lambda /\mathrm{\Delta }\lambda 400`$) of the XMM/RGS, which is not sensitive to the source extension, it is possible to overcome the above described difficulties.

M82 was observed several times by Chandra; we consider the two longest exposures, of 33 ks and 15 ks respectively, both obtained on September 1999. A cumulative spectrum of the brightest point sources was extracted from the two Chandra ACIS-I observations; the best-fit model was an absorbed power-law with $`\mathrm{N}_\mathrm{H}=7.9\pm 0.7\times 10^{21}`$ cm-2, $`\mathrm{\Gamma }=0.84\pm 0.07`$, $`\chi _\mathrm{r}^2=0.4`$, which is typical of High Mass X–ray Binaries (see Persic & Rephaeli, 2002, and reference therein). The 0.5–10 keV flux was $`4.2\times 10^{12}`$ erg s-1 cm-2, which represented $`25\%`$ of the total (point source plus diffuse) flux (but notice that M82 is a variable source, Ptak & Griffiths 1999, Rephaeli & Gruber 2002).

M82 was observed by XMM-Newton for 30 ks in May 2001; in the XMM-Newton observation, after light-curve cleaning, $`20`$ ks of data were available for scientific analysis. We extracted the EPICpn spectrum from a circular region centered on the starbursting core of M82 and with a $`15\mathrm{}`$ radius, which approximately matches the RGS PSF; the standard recipe for EPIC spectra extraction suggested in the XMM web pages was followed. The XMMSAS software version 5.3.0 was used. A background spectrum was taken from standard background files provided by the XMM-Newton SOC. RGS spectra were extracted with standard settings (i.e. including events within 90% confidence and rejecting as background events those outside 95% of the PSF figure); this procedure is similar to the one used by Read & Stevens (2002) but with more recent calibrations. After visual inspection, we discarded the 2nd-order spectra because of low signal/noise ratio. We included in the fit only the 1st-order channels comprised in the 6–20 Å interval, since outside the S/N ratio rapidly drops. A systematic error up to 10% in the 14–18 Å channels was included (XMM calibration document XMM-SOC-CAL-TN-0030 Issue 2).

The hard–band (2–7 keV) Chandra images show that the point sources emission largely exceeds that of the hot gas, while the hard part (3–7 keV) of the pn spectrum can be well described by a power-law whose best-fit slope resulted in the same value previously determined for the Chandra point sources. Thus, making use of information from both Chandra and XMM-Newton observations, we have been able to constrain the point source slope, absorption and normalization at the moment of the XMM-Newton observation.

We jointly fitted the pn and RGS spectra with XSPEC Version 11.2 (see Fig. 3) with a two component model which is the sum of: *i)* an absorbed power-law (accounting for point sources) with the slope and absorption parameters fixed at the best fit values obtained with Chandra while normalization was left free to account for variability; *ii)* an absorbed, optically–thin thermal plasma emission, with variable line intensities and differential emission measure (DEM) distribution. The DEM distribution was modeled with a 6$`^{\text{th}}`$ order Chebyshev polynomial (c6pvmkl model in XSPEC); with the best-fit parameters it resembles a bell-shaped curve with a peak at $`kT0.7`$ keV and FWHM $`0.65`$ keV (i.e. the plasma temperatures range from $``$0.35 keV to $``$1 keV), in agreement with Read & Stevens (2002). The thermal component is absorbed by cold gas with best-fit column density $`N_H3.8\times 10^{21}`$ cm-2.

The confidence contours of the Iron abundance versus normalization of the thermal component are reported in Fig. 4. As expected the two parameters are well correlated as the higher is the normalization of the continuum spectrum the lower is the intensity of the emission lines and in turn the element abundances. The derived abundances and their associated 90% confidence error are shown in Table 1. We notice that these abundances are about a factor 2 below those reported in Read & Stevens (2002); we attribute the discrepancy to the different modeling of the continuum emission (Read & Stevens 2002 did not include pn data, neither considered the contribution from point sources). However, it is worth of mentioning that both sets of abundances indicate a low O abundance compared to the other $`\alpha `$-elements.

To further check the consistency and robustness of our estimates, we repeated the fit two times considering only the RGS data and using *i)* the above described thermal plus power-law model and *ii)* a thermal model alone. This last check, in particular, most closely matches the analysis previously made by Read & Stevens (2002) to the point that the discrepancies are within the errors and can be attributed to the differences in the used instrument calibration. The abundances derived from the three different models are plotted in Fig. 5. Those obtained combining pn and RGS have significantly smaller uncertainties. One can see that different models have mainly the effect to almost rigidly scale all the abundances, with minor impact on the abundance ratios.

The low Oxygen abundance is somehow difficult to explain in the framework of $`\alpha `$-element enhancement by type II SN explosions. Thus, to check its reliability we included in the fit also the channels corresponding to the 20–23 Å interval, where the OVII triplet is observed (although with very few counts). We find no significant difference neither in Oxygen abundance nor in DEM distribution (the OVIII $`\lambda `$19/OVII $`\lambda `$22 ratio is extremely sensitive to the plasma temperature). Ne is also somewhat under-abundant. We checked for possible instrumental effects, such as variations in the RGS effective area in proximity of the Ne and O lines; however the effective area in this region is a rather smooth function of the wavelength, so that there should be no instrumental issues and these low abundances are likely to be real.

Since M82 is observed edge—on, the X-ray emission of its nuclear region is subject to quite heavy absorption which could affect the O lines detected at low energies. However, as shown in Fig. 6, the contour plot of the Oxygen abundance towards the foreground absorption, suggest a rather weak correlation. In the next section we discuss other physical possibilities to explain the O under-abundance issue.

## 5 Discussion

In M82 all the three components, namely hot and cold gas-phases and stars trace a very similar Iron abundance, the average value being \[Fe/H\]$`0.35`$ dex. Indeed, since Iron is mainly produced by SNe Ia, it is expected to be released in the ISM only after $`1`$ Gyr from the local onset of SF.

At variance, $`\alpha `$elements (O, Ne, Mg, Si, S, Ca, Ti) are predominantly released by SNe II with massive progenitors on much shorter timescales.
Stars trace an average \[$`<`$Si,Mg,Ca$`>`$/Fe\]$`0.4`$ dex, while the hot gas suggests an average \[$`<`$Si,Mg,S$`>`$/Fe\]$`0.5`$ dex. Such an overall \[$`\alpha `$/Fe\] enhancement in M82 is fully consistent with a standard chemical evolution scenario, where the ISM in the nuclear regions of massive SB galaxies is mainly enriched by the products of SNe II explosions (see e.g. Arnett 1995 for a review) occurring in recursive bursts of SF of relatively short duration.

Fig. 7 shows the various abundance ratios relative to Iron, as measured in stars, hot and cold gas-phases. Since the cold gas-phase abundances are from different data sets in the literature, hence not homogeneously determined, the inferred abundance ratios should be regarded with caution.

When comparing the hot and cold gas-phase abundances, one finds Ne over-abundant, Si and especially S significantly under-abundant in the cold gas. Förster Schreiber et al. (2001) suggests S and to a lower level Si depletion onto interstellar dust grains. However, other metals like Fe in particular, with even higher degree of incorporation into dust (Savage & Sembach, 1996), should be severely depleted in the cold gas, which does not seem the case. A possible explanation is that the optical lines used to infer Fe and O abundances (Alloin, Collin-Souffrin, & Joly, 1979) and the the mid/far IR lines used to infer Ne and S abundances (Lord et al., 2001; Förster Schreiber et al., 2001) trace different nebular sub-structures within the central few hundreds pc, somewhat chemically dishomogeneous or with different dust content.

The \[O/Fe\] abundance ratio in M82 is even more puzzling to interpret. Before the onset of the current burst of SF, the ISM was enriched in O as well as in the other $`\alpha `$-elements, as suggested by our stellar and cold gas-phase abundances. The O under-abundance measured in the hot gas cannot be modeled with a standard nucleosynthesis from the present generation of SNe II. Umeda et al. (2002) suggested explosive nucleosynthesis in core-collapse hypernovae but they predicted a large under-abundance of Ne and Mg as well, which is not observed. \[O/Fe\] under-abundance can be also explained with a major Fe enrichment by current SN-Ia explosions from previous generations of stars, but in this case all the \[$`\alpha `$/Fe\] ratios should be low, again contrary to what observed. O might be also locked into dust grains. Indeed, SNe II are known to produce significant amounts of dust (Todini & Ferrara, 2001) and O, being a basic constituent of dust grains, is expected to suffer of severe dust depletion. However, other metals like Fe in particular, but also Si and Mg should be depleted into dust (Savage & Sembach, 1996), leaving almost unchanged or eventually further enhancing the \[$`\alpha `$/Fe\] ratio.

## 6 Conclusions

Our abundance analysis of the stellar and hot gas-phase components in the nuclear region of M82 indicate an Iron abundance about half Solar and an overall $`\alpha `$enhancement by a factor between 2 and 3 with respect to Iron.

These abundance patterns can be easily explained within a standard nucleosynthesis scenario where the ISM is mainly enriched by SNe II on relative short timescales and with a star formation process occurring in recursive bursts.

Oxygen behaves in a strange fashion. It is over-abundant in stars and cold gas, similarly to the other $`\alpha `$-elements, while is significantly under-abundant in the hot gas. Major calibrations and/or modeling problems seem unlikely. Hypernovae nucleosynthesis, dust depletion, SN-Ia enrichment can somehow explain an Oxygen under-abundance but other metals should follow the same pattern, while they do not.

Somewhat exotic threshold effects and/or depletion mechanisms preferentially affecting Oxygen could be at work in the nuclear region of M82, but presently this issue remains controversial.

The financial support by the Agenzia Spaziale Italiana (ASI) is kindly acknowledged.


text:
# Implementing an Observatory Control System-I. A Generic approach

## 1 Introduction

An observatory, irrespective of wavelengths of observations, is a complex instrument consisting of a large number of interacting parts. The creation of a such a structure is a difficult problem, fraught with issues of various kinds. However, there is a great deal in common between the observatories of past, present and possibly future and it is important to create a framework which takes advantage of this commonality. In spite of the fact that scientists as a rule tend towards creating unified standardized approaches to a subject, such a standard is seemingly lacking in the astronomy related engineering and instrumentation field or is not very popular. In this paper we attempt to address the issue of consolidation of many years of experience in building such complex instruments, and creation of a framework which addresses the practical needs of the community, spanning the gamut of observatories from earth based to space based ones.

## 2 Architecture of an Observatory

The functional requirements of an observatory determines to a large extent both the physical and software architecture of the system. Unlike older observatories which had many instruments primarily interconnected via the operators, modern observatories have an internal structure which consists of parts interlinked through media capable of transmitting digital signals. The required analog conversions are done at possibly the last step, thus allowing the observatory to be modeled as a large asynchronous digital network.

Given this the problem that we are trying to address is the following: Is it possible to define a generic set of components and interconnects, starting from hardware to users of various kinds, which will allow the observatory to be modeled and built with a minimum of effort and a maximum of stable, reliable functionality. It is also desired that the structures be generic enough to permit an easy mix and match between the different observatories, while maintaining the ability to be customized to the requirements of a particular entity.

In order to address this, it is proposed that the observatory be modeled and described as a set of interacting objects (in the computer programming sense of the word) and such a model be implemented at all levels by creating hardware and software components which can be put in one to one correspondence with the model. In this paper we present one such model.

The model is intended to be an evolutionary structure, permitting instantiation of it with variations at various sites and a feedback and evolve methodology defined within.

## 3 The Common Object Framework

In order to facilitate such an architecture one needs to adopt a standard framework for designing and implementing the Observatory Control System. Moreover the standardisation attempt should extend itself to hardware (the common parts of it at least) as well as the test and calibration framework.

What is proposed here in the current paper is such a framework for a distributed observatory system. In forthcoming papers we detail the implementation aspects as applied in the IUCAA observatory project.

A modern observatory consists of interconnected digital systems from motion control to data acquisition and archiving, irrespective of wavelength or location of the system. Going further many of them are designed to operate with minimal or no manual intervention, such in the case of space observatories or robotic ones. Building and deploying such a complex system to reach level of high precision and noise free observations is currently an expensive and time consuming task. However there have been a large number of them of all sizes and categories built and deployed on Earth as well as in nearby orbits in space. The knowledge and experience that has accumulated from such successful and unsuccessful endeavours is enormous. However the knowledge has not been translated into generic patterns  which can be used to create practical and stable platforms. There is a requirement for a centralised repository of auxiliary information such as the issues faced during development and deployment as well as data regarding long term use of the instruments and so on.In the software engineering field many approaches have been developed to create patterns for distributed and concurrent systems to which class observatories certainly belong. . In a like manner we try to describe a methodology for practically constructing these systems. It has been suggested that since the technology is constantly changing and it is difficult for the astronomy community to concentrate the resources from different domains of expertise, it would be efficient to have such models implemented using currently established software development companies. It would also allow experienced resource management frameworks to be introduced into the picture for timely and reliable creation. (Prof. Ajit Kembhavi, private communication)

There are many aspects to be considered in the common framework. Other than achieving the desired functionality, the system must also be open to probing and diagnostics as well as maintenance and other management issues. Issues related to calibration and performance analysis of the system are also critical to successful deployment.

### 3.1 The basic model

The observatory is designed by assembling together various components (almost like a Lego kit). A component is defined herein as an entity consisting of state and behaviour, implemented both in hardware and software. In some senses it is a unification of hardware and software under the object oriented paradigm. The component also incorporates along with other functional aspects related to application, an interconnect architecture, including the hardware and the software protocol stack. Thus it is akin to the transputer in the world of microprocessors. In practice, the notion of “component” is defined not at the level of the controller that controls the hardware, but more at the level of the PC that interfaces to the hardware device on one side and meshes into the observatory control system via the transport layer interface. Thus an important aspect of the state of the “component” is the interconnect address. It is proposed that this communication layer be TCP–UDP/IP, the standard protocol stack of the Internet.

Over this communication link, a message based, event driven architecture is defined which allows the whole observatory to be built up by integrating these components. Self descriptive messages in XML or in alternative formats can be sent across these channels in point to point, broadcast or to a central message router and so on such that the sender and recipient can effect either synchronous or asynchronous changes of state based on them. A suite of messages relevant for this model in XML format is presented later.

Messages can have different degrees of genericity and so can the family of handlers which process the messages. There is a hierarchy of handlers starting from the most general to the most specific as in standard protocol stacks.

### 3.2 Deployment and Initialization

We begin by assuming that the final entity built out of these components have been deployed in a required spatial pattern. However we must address the issue of initialization of the same. The proposed methodology is to use a central configuration server to serve the configuration information for each component via the message protocol and update itself on receiving initialization confirmation from the components. The server is based on LDAP although any other directory/database server could be utilized for the same purpose. However the schema that is used to store and communicate the configuration information, in the form of a descriptor must be standardised. The initialization phase is modeled akin to the discovery phase of a protocol such as USB.

Initialization also involves creation of global data which can be used for communication as well as synchronization by the components. A service providing this facility for a distributed system is a component of the architecture like others.

However there is an issue of time ordering of the initialization. This is an important consideration which is addressed in the implementation stage, where the server masterminds the order in which the components are initialized inspite of the arbitrary ordering of startup in physical time.

### 3.3 The component architecture

As is clear from the above discussion, there are two degrees of freedom in design and implementation. By evolving a standard for both we should limit the variability in implementation.

We begin by designing an observatory profile which defines and characterises the components that constitute an observatory. We further define and categorise the messages that can be sent and received by the components, there by producing a constrained environment which meets our needs. The system is designed and constructed to be an open system allowing the addition of more ratified messages and components as well as leaving some free space in the object definitions and message definitions. The algorithms that the message handlers implement are also defined to a great degree by the standard but more in line with directives, rather than mandates.

## 4 The Observatory Profile

The observatory profile consists of the list of components and the messages that they transmit and receive. The profile defines the basic structure of the component and delineates an object framework for creating and classifying them. It also defines the functionality in detail along with the algorithms that are used/recommended to implement the functionality. In this section we will briefly discuss the observatory profile which is detailed in a separate paper. The profile consists of a set of other profiles which individually detail the objects and their behavioural and structural constraints.

1. Telescope System Profile: A detailed set of objects that encapsulate the telescope system and its behaviour including the autoguider and pointing models.
2. Enclosure System Profile: Set of objects that control the enclosure of the telescope system and individual instruments.
3. Instrument Profile: Set of instruments such as CCD cameras and other Imagers, Spectroscopes, Polarimeters and others.
4. Data archiving and Processing Profile: Objects to archive data in raw and processed forms (from image processing pipelines deployed on clusters and grids) such that query, retrieval etc. is uniform, distributed and simple.
5. Management System Profile: The large set of objects which are used to manage the observatory and its components such as loggers, monitors and management stations by means of which the operators can interact with the observatory system.
6. Configuration and security repository Profile: The central repository storing configuration information regarding the components of the system, the users and their authorizations and so on.
7. External Environment and Meteorology System Profile: The set of objects comprising the functionality required to access, control and monitor the external conditions such as electrical and mechanical power, weather conditions, infrastructure (networks for instance) conditions and so on.

Each individual profile is discussed in greater detail in later papers. However in the following sections we will briefly discuss some non standard subsystems.

## 5 The Management and Configuration subsystems

In order to keep such a complex system running smoothly as well as make it extensible and maintainable one must have security and management capabilities defined in the the infrastructure itself. A management model based on Common Infrastructure for Management (CIM)  over HTTP may be utilised to ensure that a stable generic management platform is integrated into the design.

In a like manner a security system suitable for such a distributed architecture must be implemented and the security credentials must be movable so that distributed resources can be accessed transparently. These issues will be discussed in detail in the later work.

## 6 The data archiving and processing system

Since the raison d’etre of the observatory is to produce data for science the data archiving and processing subsystems are extremely important. There are many different archiving structures and formats and a few established tools and patterns for astronomical data processing. One of the primary products of this track which has been evolving over time is the virtual observatory project.

The data archiving component is designed to be conformant with the ongoing efforts in this field and will be detailed in a later paper. The retrieval, processing and caching of data required for creating distributed, streaming mode architectures suitable for parallel or grid based data processing and integrating the same into existing applications such as IRAF and MIDAS will also be discussed later.

## 7 Conclusions

In this paper we present the principles and architecture of a distributed observatory control system in a generic form. We propose that such a standardised componentised architecture be used to simulate, prototype and build the observatories in all wavelengths. We propose a language independent implementation of the objects in the object oriented framework which describes both the components as well as the gluing structures as objects. In succeeding papers we discuss the actual implementation in greater detail.

## 8 Acknowledgements

The author wishes to thank IUCAA and his colleagues in the instrumentation division for much help in the implementation of this architecture. Many issues and problems were understood and resolved during discussions with Prof. S.N. Tandon of IUCAA. We also benefited greatly from extended discussions on many aspects with Prof. Padmanabhan of IUCAA and we thank him for all the guidance.


text:
# Comment on “Accurate ground-state phase diagram of the one-dimensional extended Hubbard model at half filling”

## Abstract

It is shown that Guoping Zhang’s results \[G.P. Zhang, Phys. Rev. B68, 153101 (2003)\] for the charge-density-wave phase boundary in the half-filled one-dimensional extended Hubbard model are incorrect and that his criticism of my work \[E. Jeckelmann, Phys. Rev. Lett. 89, 236401 (2002)\] is groundless.

In Ref. Zhang, , Guoping Zhang presented density-matrix renormalization group (DMRG) results which contradict my DMRG calculations Jeckelmann  and Hirsch’s quantum Monte Carlo (QMC) simulations Hirsch  for the charge-density-wave (CDW) phase boundary in the one-dimensional extended Hubbard model at half filling. In this Comment I show that Guoping Zhang’s results are inaccurate and that his criticism of my work is groundless.

Although the phase diagram of the extended Hubbard model is still partially controversial (see Refs. Jeckelmann, ; Sandvik, ; YZZhang,  and references therein), the CDW phase boundary $`V_c(U)`$ in the parameter space ($`U,V`$) was determined years ago Hirsch ; Cannon  and has not been disputed in recent studies. Jeckelmann ; Sandvik ; YZZhang ; Sengupta ; Nakamura  In Fig. 1, I show the results of various numerical investigations for $`V_c(U)U/2`$ in the weak to intermediate coupling regime. There is an excellent overall agreement between Hirsch’s QMC simulations, Hirsch  the exact diagonalizations of Cannon et al.Cannon  Nakamura’s level crossing analysis, Nakamura  the stochastic series expansion QMC (SSE-QMC) simulations of Sandvik et al.Sengupta ; Sandvik  Yuzhong Zhang’s DMRG calculations YZZhang , and my DMRG calculations Jeckelmann . In particular, my results agree quantitatively errors  with the most recent and accurate numerical simulations. Nakamura ; Sengupta ; Sandvik ; YZZhang  Only Guoping Zhang’s DMRG data Zhang ; oldZhang  deviate systematically from the other results. Therefore, there is clearly a problem with his calculations.

The discrepancy between the various DMRG calculations Zhang ; Jeckelmann ; YZZhang ; oldZhang  is not surprising. Guoping Zhang uses the infinite-system DMRG algorithm while Yuzhong Zhang and I use the more accurate finite-system DMRG algorithm. White  It is well known Schollwoeck  that for many problems the infinite-system algorithm yield incorrect results while the finite-system algorithm gives essentially exact (numerical) results. In particular, it is essential to use the more reliable finite-system DMRG algorithm for inhomogeneous systems such as a CDW ground state. Therefore, the discrepancy between Guoping Zhang’s results and all other works just demonstrates the failure of the standard infinite-system DMRG algorithm for the present problem. (See Ref. Prokofiev,  for another example of the infinite-system algorithm failure and Ref. Rapsch,  for the successful investigation of the same problem with the finite-system algorithm.)

In his paper Guoping Zhang wrongly claimed that my DMRG calculations (and the QMC simulations of Ref. Hirsch, ) failed to reproduce the weak coupling limit result $`V_c(U)=U/2`$. In Refs. Jeckelmann, ; Hirsch,  the investigation of the phase diagram was focused on the intermediate- and strong-coupling regimes (i.e., $`U2t`$) and no analysis of the weak-coupling limit $`Ut`$ was performed. Here I present additional results for $`V_c(U)`$ calculated with DMRG for weaker couplings: $`V_c/t=0.260\pm 0.003`$ for $`U/t=0.5`$, $`V_c/t=0.545\pm 0.005`$ for $`U/t=1`$, and $`V_c/t=0.835\pm 0.005`$ for $`U=1.5t`$. Moreover, I have calculated $`V_c(U=2t)/t`$ more accurately and found $`1.113\pm 0.005`$ (in agreement within the error bars with the value given in Ref. Jeckelmann, ). These results are shown in Fig. 1 as square. They agree perfectly with other works Sandvik ; YZZhang ; Nakamura  and, clearly, they approach the weak-coupling result $`V_c(U)=U/2`$ in the limit $`U0`$. In Figs. 1 and 2 of Ref. Zhang,  Guoping Zhang used a different representation of the data, $`U/V_c(U)`$ vs $`U`$, to analyze the weak-coupling limit. In Fig. 2, I show again all data of Fig. 1 using this representation. Clearly, the minimum of the ratio $`U/V_c(U)`$ occurs for $`U`$ slightly smaller than $`2t`$ and the weak-coupling limit is recovered only for $`U`$ smaller than $`2t`$. Therefore, the $`U0`$ limit of $`U/V_c(U)`$ cannot be determined using numerical data for $`U2t`$ and the Figs. 1 and 2 of Ref. Zhang,  are misleading. The DMRG and QMC data for $`U2t`$ presented in Refs. Jeckelmann,  and Hirsch,  are fully compatible with the weak-coupling limit $`V_c(U)=U/2`$, contrary to Guoping Zhang’s assertion in Ref. Zhang, .

In summary, comparisons with the results available in the literature confirm the accuracy and reliability of the DMRG calculations presented in Ref. Jeckelmann, . Guoping Zhang’s results and conclusion are faulty due to the inappropriate use of the infinite-system DMRG algorithm.


text:
# Gravitomagnetic measurement of the angular momentum of celestial bodies

## Abstract

The asymmetry in the time delay for light rays propagating on opposite sides of a spinning body is analyzed. A frequency shift in the perceived signals is found. A practical procedure is proposed for evidencing the asymmetry, allowing for a measurement of the specific angular momentum of the rotating mass. Orders of magnitude are considered and discussed.

A well known effect of gravity on the propagation of electromagnetic signals is the time delay: for example an electromagnetic beam emitted from a source on Earth toward another planet of the solar system, and hence reflected back, undergoes a time delay during its trip (with respect to the propagation in flat space-time), due to the influence of the gravitational field of the Sun. This effect was indeed called the ”fourth” test of General Relativity,, coming fourth after the three classical ones predicted by Einstein.

Though small, the time delay in the propagation of electromagnetic waves was detected by Shapiro et al. , timing radar echoes from Mercury and Venus, by means of the radio-telescopes of Arecibo and Haystack. Anderson et al.  measured the time delay of the signals transmitted by Mariner 6 and 7 orbiting around the Sun. Finally Shapiro and Reasenberg obtained more accurate results using a Viking mission that deposited a transponder on the surface of Mars: the theoretical prediction was verified within $`\pm 0.1\%`$  .

These measurements accounted just for the presence of a massive source, described by the Schwarzschild solution of the Einstein field equations. However, there is another correction to the time delay, due to the spin of the source, which is produced by the so called ”gravitomagnetic” interaction.

The time delay in the gravitational field of the Sun can be studied in the weak field approximation. By using Cartesian coordinates and assuming that the $`z`$ axis coincides with the direction of the angular momentum of the source, the propagation of electromagnetic signals, in the equatorial plane, is described by the null world-line

$$0=g_{tt}dt^2+g_{xx}dx^2+g_{yy}dy^2+2g_{xt}dxdt+2g_{yt}dydt$$
(1)
In the weak field approximation, the bending of the trajectories due to the gravitational field is neglected. This assumption, appropriately choosing the $`x`$ and $`y`$ axes, leads to a ray trajectory that is a straight line $`x=b=const`$; $`b`$ is of course the closest approach distance with respect to the spinning body. On using the appropriate form of the metric elements, the time of flight of the electromagnetic beams can be written as

$$t_f(y_1,y_2)=t_0+t_M+t_J$$
(2)
where

$`t_0`$ $`=`$ $`{\displaystyle \frac{y_2y_1}{c}}`$ (3)
$`t_M`$ $`=`$ $`{\displaystyle \frac{2GM}{c^3}}\mathrm{ln}{\displaystyle \frac{y_2+\sqrt{b^2+y_2^2}}{y_1+\sqrt{b^2+y_1^2}}}`$ (4)
$`t_J`$ $`=`$ $`{\displaystyle \frac{2GJ}{c^4b}}\left[{\displaystyle \frac{y_2}{\sqrt{b^2+y_2^2}}}{\displaystyle \frac{y_1}{\sqrt{b^2+y_1^2}}}\right]`$ (5)
In Eqs. (4) and (5) $`M`$ is the mass of the source and $`J=Mca`$ is its angular momentum, assumed to be orthogonal to the plane of the motion.

The quantity $`y_1`$ is the $`y`$ coordinate of the source of the signals, $`y_2`$ is the $`y`$ coordinate of the receiver. The time $`t_0`$ is clearly the Newtonian time of flight; $`t_M`$ is the gravitational time delay measured by Shapiro et al., and $`t_J`$ is the correction to the time delay produced by the gravitomagnetic interaction with the angular momentum of the central body. The double sign in $`t_J`$ means that gravitomagnetism shortens the time of flight on the left and lengthens it on the right.

This asymmetry, within the solar system, is in any case small: however its systematicness lends an opportunity to reveal the effect, appropriately combining the ticks of a ’clock’ passing behind the spinning mass. Elsewhere we have outlined a way to reveal this effect, based on the fact that the relative motion of source, receiver and central mass, produces a varying time delay, which shows up as a small frequency shift. The gravitomagnetic contribution to this shift is manifested as an asymmetry between right and left with respect to the central mass. Let $`\nu _0`$ be the proper frequency of the electromagnetic signal; in these conditions, mirroring the signal after the occultation (time reflection), then superposing corresponding records with the signal before the occultation of the clock, a beating function will result, where the magnitude of the frequency of the beats is proportional to the angular momentum of the spinning body

$$\nu _2=4\nu _0\frac{\mu a}{b^2}\frac{v_0}{c}$$
(6)
and the frequency of the basic signal is shifted with respect to the flat space-time situation by an amount proportional to the mass:

$$\nu _1=\nu _0\left(1+4\frac{\mu }{b}\frac{v_0}{c}\frac{\mu b}{R^2}\frac{v_0}{c}\right)$$
(7)
In Eqs. (6) and (7) $`v_0`$ is the apparent transverse velocity of the source in the sky, $`R`$ is the distance of the observer from the central body and $`\mu =G\frac{M}{c^2}`$. To fix a few numbers, let us consider the situation in the solar system, with the Sun as the central spinning body, and an Earth bound observer. The source is a far away astronomical body (e.g. a pulsar). In this case the orders of magnitude are

$$\begin{array}{ccccc}\mu 10^3m& a10^3m& b10^9m& R10^{11}m& v_010^4m/s\end{array}$$

As a consequence we have

$`4{\displaystyle \frac{\mu }{b}}{\displaystyle \frac{v_0}{c}}`$ $``$ $`10^{10}`$
$`{\displaystyle \frac{\mu b}{R^2}}{\displaystyle \frac{v_0}{c}}`$ $``$ $`10^{14}`$ (8)
$`4{\displaystyle \frac{\mu a}{b^2}}{\displaystyle \frac{v_0}{c}}`$ $``$ $`10^{16}`$
The frequency effect, connected with a varying time delay in the propagation of electromagnetic signals in the vicinity of a spinning massive body, is very small, however not entirely negligible, at least when a suitable procedure to detect it is used, such as the one we have just outlined, based on the possibility of producing a beat between signals after and before the occultation.

Of course many practical problems have to be considered and discussed to transform some principle formulae into an actual measurement. Other situations of physical interest can be taken into account: for instance, a more favorable observational condition (with regard to the amount of the effect) can be the one of a pulsar orbiting a neutron star. For more details, see the discussion in our paper. What we can say, at this stage, is that the method we propose seems promising for detecting gravitomagnetic effects on the propagation of electromagnetic signals, and for measuring the angular momenta of celestial bodies.


text:
# Spin-boson models for quantum decoherence of electronic excitations of biomolecules and quantum dots in a solvent

## 1 Introduction

Biomolecular systems contain a diverse range of optically active molecules that are crucial to biological function . Important examples include retinal in rhodopsin which plays a role in vision, porphyrins and chlorophyls which play roles in photosynthesis , photoactive yellow protein, blue copper proteins, and the green fluorescent protein (GFP)  which has become a powerful tool in molecular and cell biology. The different chromophores couple very differently to their environment which consists of the surrounding proteins and the solvent (water or aqueous electrolytes) leading to very different quantum dynamics . Related questions arise concerning the spatial extent of the quantum coherence of excitations in spatially extended structures such as $`\alpha `$-helices and $`\beta `$-sheets in proteins and rings of porphyrin molecules in light harvesting complexes . Thus, it is important to understand the coupling of electronic excitations in chromophores to the environment. This is not just of interest for understanding biological function, but also because these biomolecular systems could potentially be model systems to understand quantum decoherence at the nanoscale. This is because they have a complexity that cannot be fabricated even with the most advanced nanotechnology, but they can be modified and tuned using sophisticated biochemical techniques . We also note that there are potential applications in biosensors with biomolecules coupling to quantum dots via resonant energy transfer , and with biomolecules acting as optically controllable elements in molecular electronics  or as qubits in a quantum computer .

Understanding the quantum dynamics of a system that is strongly coupled to its environment is a highly non-trivial, and sometimes controversial, problem. Furthermore, given the complexity and diversity of biomolecular systems it is desirable to find tractable model Hamiltonians which capture the essential physics and are amenable to rigorous analysis. The past few decades have seen monumental advances in combined quantum mechanical and molecular dynamics simulations, allowing the study of biomolecular systems of realistic size. See for example, the recent study of the quantum dynamics of retinal in the rhodopsin environment . However, in spite of this progress it is unrealistic that in the next decade it will be possible to give a full treatment of the quantum dynamics including the quantum many-body effects associated with the strong interaction between the chromophore and its environment. Furthermore, simple models amenable analytical treatment can give insight into the essential physics involved with qualitative differences, especially when model parameters are extracted from quantum chemistry and molecular dynamics.”

We first consider the problem of the interaction of an electronic excitation (exciton) in a biomolecule (or a quantum dot) with the surrounding polar solvent. We look for the simplest possible model which can still describe the relevant physics, which allows us to determine the important parameters and perform analytical calculations. We show that this system can be described by an independent boson model where the environment is modelled as an infinite collection of harmonic oscillators with ohmic response. In the model we derive, the strength of the coupling to the environment is determined by the frequency dependent dielectric function of the polar solvent and by the difference between the electric dipole moments of the molecule in the ground and excited states. The relaxation rate of the polar molecules within the solvent determines the cutoff frequency for the ohmic response.

We then introduce minimal models which can describe the effect of the solvent on the quantum dynamics associated with conformational changes at conical intersections , and on excitations which are coupled via a mechanism such as resonance energy transfer. The latter is important in photosynthesis , is the basis of Fluoresence Resonance Energy Transfer (FRET) spectroscopy which is used to determine distances between chromophores in biomolecules , and is the basis of new biosensors . We have shown that each may be modelled by spin-boson models  which exhibit rich many-body physics, and have been used to model systems such as the coupling between electron transfer and protein motion or a solvent . The general spin-boson Hamiltonian is

$$H=\frac{1}{2}ϵ\sigma _z\frac{1}{2}\mathrm{\Delta }\sigma _x+\underset{\alpha }{}\omega _\alpha a_\alpha ^{}a_\alpha +\frac{1}{2}\sigma _z\underset{\alpha }{}C_\alpha (a_\alpha +a_\alpha ^{}),$$
(1)
where $`\sigma _z,\sigma _x`$ are Pauli matrices describing the TLS, $`ϵ`$ is the energy separation of the two levels, $`\mathrm{\Delta }`$ is the tunneling matrix element between the two levels, and the bath is modelled by boson operators $`a_\alpha ,a_\alpha ^{}`$ with frequencies $`\omega _\alpha `$ and couplings $`C_\alpha `$ to the TLS. The independent boson model is (1) with $`\mathrm{\Delta }=0`$. This means that, unlike the spin-boson model, there are no transitions between the two quantum states. Only one environmental parameter, however, contributes to the system dynamics, the *spectral density*

$$J(\omega )=4\pi \underset{\alpha }{}C_\alpha ^2\delta (\omega \omega _\alpha ).$$
(2)
Leggett *et al*.  defined a dimensionless coupling constant $`\alpha `$ which determines the extent of the quantum coherence. We show that for realistic biomolecules solvated in water this parameter will typically have values of order unity, which would prevent coherent Bloch oscillations for degenerate electronic levels. Such large values are in stark contrast to the Josephson junction qubits for which the effect of the Johnson noise in the electronic circuit is described by a spin boson model with values of $`\alpha `$ that are typically many orders of magnitude smaller .

## 2 Independent boson model for the interaction of a single chromophore with a solvent.

We now show how an independent boson model can describe an electronic excitation in a biomolecule or quantum dot coupled to a solvent bath. We model the chromophore as a TLS with only a single active transition, of energy $`ϵ`$. We assume the solute molecule has a permanent dipole moment in both the ground and excited state, given by $`\stackrel{}{\mu }_g`$ and $`\stackrel{}{\mu }_e`$ respectively. Describing the electronic ground and excited states by fermion creation operators $`c_g^{}`$ and $`c_e^{}`$ respectively, the solute part of the total Hamiltonian operator is

$$H_{TLS}=\frac{1}{2}ϵ(c_e^{}c_ec_g^{}c_g).$$
(3)
In the simplest picture, we can describe the solvation process by the well-known Onsager model . This is a type of *continuum model* where the solvent is treated as a homogeneous dielectric, and has been used for both analytic and computational studies of solvation . The solute is treated as a point dipole which is surrounded by a spherical cage of polar solvent molecules with Onsager radius $`a`$ (Figure 1), which is typically the size of the solute molecule . The cavity is assumed to be a vacuum, i.e., it has a dielectric constant $`ϵ_r=1`$. The hydration shell of water molecules who have broken hydrogen bonds with their neighbours, and are now bonded to the biomolecule, will serve to increase the effective cavity radius $`a`$. Further details (e.g., surrounding proteins) can be treated by adding shells with different dielectric constants around this cavity; this will be considered in a later work.

The central dipole polarises the cage, which in turn produces an electric field inside the cavity, called the *reaction field* $`\stackrel{}{R}(t)`$. For the case of a uniform spherical cavity the reaction field is constant everywhere inside . This field acts back on the dipole, with interaction energy $`E=\stackrel{}{\mu (t)}\stackrel{}{R(t)}`$, typically lowering the total energy and hence forming a stable structure. Although the Onsager model does not include the microscopic details of the system, it does capture the essential physics of the solvation process  and provided due care is taken to be aware of important microscopic effects (such as charge transfer processes) it should still provide valuable physical insight.

Prior to the electronic transition, we assume the system to be in thermal equilibrium, and the reaction field will be parallel to the dipole $`\stackrel{}{\mu }_g`$. For simplicity, we assume that the direction of the dipole moment of the solute molecule in the ground and excited states points in the same direction. Defining $`\mathrm{\Delta }\mu =\mu _e\mu _g`$, $`\sigma _zc_e^{}c_ec_g^{}c_g`$, and identifying $`c_e^{}c_e+c_g^{}c_g`$ with the $`2\times 2`$ identity,

$`H_{solute}={\displaystyle \frac{1}{2}}ϵ\sigma _z,`$ $`H_{int}`$ $`={\displaystyle \frac{1}{2}}\sigma _z(\mathrm{\Delta }\mu )R+{\displaystyle \frac{1}{2}}(\mu _e+\mu _g)R.`$
For a quantum mechanical model, we must quantize the reaction field. In the Heisenberg picture, we express the reaction field $`R(t)`$ in its Fourier transform modes,

$$R(t)=\underset{\alpha }{}e_\alpha \left[a_\alpha e^{i\omega _\alpha t}+a_\alpha ^{}e^{i\omega _\alpha t}\right],$$
(4)
and quantise the coefficients so that $`a_\alpha `$ and $`a_\alpha ^{}`$ obey the boson commutation relations $`[a_\alpha ,a_\beta ^{}]=\delta _{\alpha ,\beta }`$.

We are therefore modelling the environment (solvent) as a bath of independent harmonic oscillators. The energy stored in the solvent cage is now expressed as $`H_{solvent}=_\alpha \omega a_\alpha ^{}a_\alpha `$ and the full Hamiltonian for a single biomolecule coupled to the environment can then be written as

$`H`$ $`={\displaystyle \frac{1}{2}}ϵ\sigma _z+{\displaystyle \underset{\alpha }{}}\omega _\alpha a_\alpha ^{}a_\alpha +`$ (5)
$`\sigma _z{\displaystyle \underset{\alpha }{}}M_\alpha (a_\alpha ^{}+a_\alpha )+{\displaystyle \underset{\alpha }{}}\stackrel{~}{M}_\alpha (a_\alpha ^{}+a_\alpha ),`$
where we define couplings $`M_\alpha =\frac{1}{2}e_\alpha \mathrm{\Delta }\mu `$ and $`\stackrel{~}{M}_\alpha =\frac{1}{2}e_\alpha (\mu _e+\mu _g)`$.

This is the independent boson model . We have not included a $`\sigma _x`$ term, so that there can be no transitions between the two energy levels of the chromophore and the chromophore energy $`ϵ\sigma _z`$ is a constant. Two conditions are needed to justify this. First, we consider only times after the initial excitation of the system (e.g., by a laser pulse) but much shorter than the radiative lifetime of the chromophore (typically 1-10ns ). This means we do not include coupling of the TLS to the quantum and thermal fluctuations of the electric field which can induce optical transitions. At energies of the order of eV these fluctuations are much smaller than the dielectic fluctuations associated with the solvent. Secondly, for many electronic excitations we have assumed that the energy gap $`ϵ`$ (typically of order of 1eV) is significantly larger than the tunneling matrix element, and so its effects can be ignored on the timescales of interest. We defer discussion of it until below when we consider $`ϵ`$ becoming small near a conformational transition.

The independent boson model can be used to investigate decoherence of systems coupled to a bath. Consider an initially coherent excitation, i.e., a coherent superposition of the ground and excited states of the chromophore. In the reduced density matrix of the chromophore, the population of each level is given by the diagonal elements, while the coherence is given by the off-diagonal elements. Thus the initial reduced density matrix of the chromophore has non-zero off-diagonal terms. Reference  shows that while the populations remain constant (as expected) the coherence decays exponentially, as given by Equations (18) and (19) and Figure 7 of . Hence, after long times the reduced density matrix has no diagonal terms and so describes the chromophore in a mixed state .

The above model therefore describes the decoherence of an initially coherent excitation of a chromophore due to its environment (water, proteins), as well as spectral features such as the shifting and broadening of the absorption peak, as will be described below.

## 3 Derivation of the spectral density function

To complete the model we must specify the environmental coupling, i.e., the spectral function $`J(\omega )`$ from (2). Following an approach similar to Caldeira and Leggett  who considered a Josephson junction coupled to an electronic circuit, we relate $`J(\omega )`$ to the zero-temperature fluctuations in the uncoupled environment . Noting that with no solute-solvent interaction $`R(t)=0`$, we now examine the reaction field fluctuation correlation function $`S(t)`$  defined as

$$S(t)=iR(t)R(0)\theta (t)i0\left|e^{iHt}Re^{iHt}R\right|0\theta (t),$$
(6)
where $`|0`$ is the ground state of the solvent harmonic oscillators, and $`HH_{solvent}=_\alpha \omega _\alpha a_\alpha ^{}a_\alpha `$. We shall see that it is the imaginary part of the Fourier transform of $`S(t)`$, which we denote $`(\omega )`$, that is relevant. We can write this as a sum over energy eigenstates ,

$$\mathrm{Im}(\omega )=\pi \underset{n}{}\delta (\omega E_n)\left|0\left|R\right|n\right|^2.$$
(7)
Expanding $`R`$ into its normal modes again, we see that all terms $`0\left|R\right|n`$ vanish except for when a single oscillator is singly occupied. These states have energy $`E_n=\omega _\alpha `$, and $`0\left|R\right|\alpha =e_\alpha `$. Therefore, we see $`\mathrm{Im}(\omega )=\pi _\alpha e_\alpha ^2\delta (\omega \omega _\alpha )`$, and using (2), we find

$$J(\omega )=(\mathrm{\Delta }\mu )^2\mathrm{Im}(\omega ).$$
(8)
To calculate the reaction field fluctuations, we note that in the Onsager model , $`R(t)`$ and the central dipole $`\mu (t)`$ are related by a linear response function $`\chi (tt^{})`$ , such that in Fourier space $`R(\omega )=\chi (\omega )\mu (\omega )`$. The static limit, where the solvent cage adjusts instantaneously, so that $`R(t)=\chi _s\mu (t)`$, is well known . In reality, the solvent cage will lag behind a changing dipole due to dielectric friction ,and $`\chi (\omega )`$ is given  by

$$\chi (\omega )=\frac{1}{4\pi ϵ_0a^3}\frac{2(ϵ(\omega )1)}{2ϵ(\omega )+1}.$$
(9)
where $`ϵ(\omega )`$ is the frequency dependent dielectric constant of the solvent. We use the Debye formula

$$ϵ(\omega )=ϵ_{\mathrm{}}+\frac{ϵ_sϵ_{\mathrm{}}}{1i\omega \tau _D},$$
(10)
where $`ϵ_{\mathrm{}}`$ and $`ϵ_s`$ are the high and low (static) frequency limits respectively, and $`\tau _D`$ is the Debye relaxation time, the bulk reorientational relaxation time of the solvent dipoles . For water, these are given by $`ϵ_s=78.3`$, $`ϵ_{\mathrm{}}=4.21`$ and $`\tau _D=8.2\mathrm{ps}`$ .

We then apply the quantum fluctuation-dissipation relation , for the imaginary part of $`\chi (\omega )`$

$$\chi ^{\prime \prime }(\omega )=(i/2)\left(1e^{\beta \omega }\right)(\omega ),$$
(11)
which reduces at zero temperature to $`(\omega )=2i\chi ^{\prime \prime }(\omega )`$. Note that this limit is well defined, even for $`\omega =0`$. Using (10) and (8), the spectral density is then given by

$$J(\omega )=\frac{(\mathrm{\Delta }\mu )^2}{2\pi ϵ_0a^3}\frac{6(ϵ_sϵ_{\mathrm{}})}{(2ϵ_s+1)(2ϵ_{\mathrm{}}+1)}\frac{\omega \tau _E}{\omega ^2\tau _E^2+1},$$
(12)
where $`\tau _E=\frac{2ϵ_{\mathrm{}}+1}{2ϵ_s+1}\tau _D`$.

This form for $`J(\omega )`$, which has a specific microscopic basis, is an important new result and could be used as an input into more phenomenological models such as the quantum Brownian model . In specific situations the magnitude of the $`J(\omega )`$ introduced here should also be compared to the $`J(\omega )`$ due to slow protein motion, as introduced in  and .

We have thus shown that solute-solvent interaction in the Onsager picture can be modeled by a independent boson model with spectral density (12). We note that $`J(\omega )`$ has an implicit high frequency cut-off at $`\omega _c=\frac{1}{\tau _E}`$, related to the finite relaxation time of the solvent dipoles. Below the cut-off, $`J(\omega )`$ is approximately linear: $`J(\omega )=\eta \omega `$, where

$$\eta =\frac{(\mathrm{\Delta }\mu )^2}{4\pi ϵ_0a^3}\frac{6(ϵ_sϵ_{\mathrm{}})}{(2ϵ_s+1)^2}\tau _D.$$
(13)
Thus, our model falls into the class of models defined in Ref.  which below some high frequency cutoff have Ohmic dissipation. $`\eta `$ should correspond to the classical frictional coefficient. Indeed, $`\eta `$ is identical to the Nee-Zwanzig prediction of dielectric friction .

## 4 Derivation of the optical absorption lineshape

We can also use this model to derive an expression for the absorption lineshape $`\alpha (\omega )`$ of the solute molecule in the presence of the solvent. In the gas phase, i.e., without the solvent, we would expect a sharp peak in the absorption spectrum at $`\omega =ϵ`$. In the presence of the solvent, however, this peak is shifted and broadened  by the solute-solvent dipole-dipole interactions.

We use the results obtained for the *independent boson model* . The resulting spectrum is

$`\alpha (\omega )`$ $`={\displaystyle \frac{A}{\omega }}\mathrm{Re}{\displaystyle _{\mathrm{}}^{\mathrm{}}}\mathrm{dt}\mathrm{exp}\left\{\mathrm{it}\left[\omega ϵ+\stackrel{~}{\mathrm{\Delta }}\mathrm{g}(\mathrm{t})\right]\right\},`$ (14)
where $`A`$ is a normalisation factor (which we are not concerned with), and $`\stackrel{~}{\mathrm{\Delta }}=4_\alpha \frac{M_\alpha \stackrel{~}{M}_\alpha }{\omega _\alpha }`$. $`g(t)`$ is the factor responsible for broadening the lineshape . We express it  using the spectral density (2):

$$g(t)=𝑑\omega \frac{J(\omega )}{\pi \omega ^2}\left[\mathrm{coth}\left(\frac{\beta \omega }{2}\right)\left(1\mathrm{cos}\left(\omega t\right)\right)+i\mathrm{sin}\left(\omega t\right)\right].$$

Since $`J(\omega )0`$ for $`\omega \omega _c`$, we can restrict the integral to $`\omega <\omega _c`$, and note that for water at room temperature (300 K, $`\tau _D8`$ ps ) $`\beta \omega <\beta \omega _c1`$. We also make the short-time approximation  and Taylor expand the trigonometric functions to give $`g(t)=\frac{1}{2}k_BT\lambda _0t^2+\frac{i}{2}\lambda _0t`$, where we have defined the reorganisation energy $`\lambda _0=\frac{2}{\pi }\frac{J(\omega )}{\omega }𝑑\omega `$. We note also that because $`\stackrel{~}{M}_\alpha =\frac{\mu _e+\mu _g}{\mathrm{\Delta }\mu }M`$, we also have $`\stackrel{~}{\mathrm{\Delta }}=\frac{\mu _e+\mu _g}{2(\mathrm{\Delta }\mu )}\lambda _0`$. Performing the Gaussian integral, we obtain the absorption spectrum

$$\alpha (\omega )=\frac{A}{\omega \sqrt{k_BT\lambda _0}}\mathrm{exp}\left\{\left[\frac{\omega \left(ϵ(\mu _g/\mathrm{\Delta }\mu )\lambda _0\right)}{\sqrt{2k_BT\lambda _0}}\right]^2\right\},$$

where we have renormalised $`A`$, and

$$\lambda _0=\frac{(\mathrm{\Delta }\mu )^2}{4\pi ϵ_0a^3}\frac{6(ϵ_sϵ_{\mathrm{}})}{(2ϵ_s+1)(2ϵ_{\mathrm{}}+1)}.$$
(15)
If we set $`ϵ_{\mathrm{}}=1`$ then the resulting lineshape is precisely that of , where the absorption spectrum was obtained classically by assuming a Boltzmann distribution of the environment. Thus, our quantum model gives the correct classical limit.

## 5 Spin-boson models of conformational change through conical intersections

This model of solvent interactions can be extended to other situations. A number of important photochemical reactions in biomolecules involve a conformational change that initiates a signal. Examples include the *cis*/*trans* isomerisation of retinal  and the photoactive yellow protein (PYP) . The former drives the proton pumping (in rhodopsin) against an external pH gradient that leads to the signal responsible for vision. The latter causes a conformational transition of the surrounding PYP protein. Only recently has it been shown  that these isomerisations occur through a “conical intersection” , between the potential energy surfaces of the ground and excited states.

A minimal model for this requires two nuclear co-ordinates (degrees of freedom) , which form a two-dimensional simple harmonic oscillator, corresponding to the classical “cone shaped” potential surfaces . The simplest model, with classical nuclear co-ordinates, which can describe a conical intersection is

$$H=ϵ\sigma _z+\mathrm{\Delta }\sigma _x+\frac{1}{2m_a}(p_a^2+m_a^2\omega _a^2q_a^2)+\frac{1}{2m_b}(p_b^2+m_b^2\omega _b^2q_b^2)+g_a\sigma _zq_a+g_b\sigma _xq_b,$$
(16)
where $`q_{a,b}`$ and $`p_{a,b}`$ represent the position and momentum respectively of the two classical oscillators, $`\omega _{a,b}`$ and $`m_{a,b}`$ are the frequencies and masses, respectively, of the oscillators, and $`g_{a,b}`$ is the coupling of each oscillator to the chromophore. This model includes a coupling term $`\mathrm{\Delta }\sigma _x`$ between the two states, where $`\mathrm{\Delta }`$ is the matrix element responsible for the non-adiabatic mixing of the ground and excited states . The spectrum is

$$E_\pm =\pm \sqrt{\left(ϵ+g_aq_a\right)^2+\left(\mathrm{\Delta }+g_bq_b\right)^2}+\frac{1}{2}(\omega _a^2p_a^2+q_a^2+\omega _b^2p_b^2+q_b^2)$$
(17)
which clearly has a conical intersection when the first term vanishes (Figure 2). We note that for $`ϵ=0`$ this corresponds to the $`Eϵ`$ Jahn-Teller model . Hahn and Stock  recently studied such a minimal model for retinal in rhodopsin but did not include the effect of the solvent, which can play a pronounced role on the topology of the conical intersection . The generalisation of this model to include the solvent and quantum motion of the nuclear coordinates is

$`H={\displaystyle \frac{1}{2}}ϵ\sigma _z+\mathrm{\Delta }\sigma _x+\sigma _z{\displaystyle \underset{\alpha }{}}M_\alpha (a_\alpha ^{}+a_\alpha )+{\displaystyle \underset{\alpha }{}}\omega _\alpha a_\alpha ^{}a_\alpha +`$ (18)
$`\omega _aa^{}a+\omega _bb^{}b+g_a\sqrt{{\displaystyle \frac{2\mathrm{}}{m_a\omega _a}}}\sigma _z(a+a^{})+\sqrt{{\displaystyle \frac{2\mathrm{}}{m_b\omega _b}}}g_b\sigma _x(b+b^{}),`$ (19)
where $`a^{}`$ and $`b^{}`$ are creation operators for the vibrations for the two degrees of freedom associated with the conical intersection.

An important question about the quantum dynamics of the chromophore retinal concerns the explanation for the following observation. The speed and quantum efficiency of the conformational change is much greater when the retinal is in the rhodopsin environment rather than in a solvent. This model could be used to address this question by performing the following (non-trivial) calculation. Suppose following optical excitation the TLS is in the excited state and the bath and vibrational modes are in the equilibrium state associated with the ground state of the TLS. We want to calculate the branching ratio between the two possible final states of the system: to the ground electronic state and the excited electronic state and the associated phonon equilibrium state. Transition to the former requires that the non-adiabatic coupling term induce a transition to the ground state by a mechanism similar to Landau-Zener tunneling. However, in the presence of strong decoherence this tunneling will be suppressed.

We now consider the magnitude and effect of the coupling of the bath to the TLS near the conical intersection. In the case of Ohmic dissipation, we can define a dimensionless coupling $`\alpha =\eta /2\pi \mathrm{}`$. For the spin-boson model (1) with $`\mathrm{\Delta }\mathrm{}\omega _c`$, if $`\alpha >\frac{1}{2}`$ and $`ϵ=0`$ there will always be incoherent relaxation , while for $`\alpha >1`$ and small $`\mathrm{\Delta }`$ the system is localised at its initial state. For $`ϵ=0`$ and $`\alpha <\frac{1}{2}`$, coherent oscillations between the two states will occur. These will be inhibited by increasing the bias $`ϵ`$ and the temperature  but enhanced by values of $`\mathrm{\Delta }`$ much larger than $`\omega _c`$ .

We can make an estimate for $`\alpha `$ in water, using the constants given earlier which completely describe the solvent. The two solute unknowns are the cavity size and change in dipole moment. If we measure these in angstroms and Debye respectively ($`1D=3.3\times 10^{30}Cm`$) then $`\alpha =22\frac{(\mathrm{\Delta }\mu )^2}{a^3}`$. For a typical chromophore with radius between $`58\text{Å}`$, a dipole moment change of just $`15D`$ is sufficient to make $`\alpha >1`$, and this condition seems likely to be met for most small molecules \[56, Ch. 6\]. In contrast, we note that for GFP, quantum chemistry calculations suggest $`\mathrm{\Delta }\mu `$ is small . The solvent is therefore a large potential source of decoherence which destroys any quantum coherent oscillations of the electronic excitation.

## 6 Fluorescence resonant energy transfer

Another important process in biomolecules is coherent transfer of excitons. An example is fluorescence resonance energy transfer (FRET), as shown in Figure 3. This occurs in photosynthesis  and is the basis of FRET spectroscopy  and a proposal for quantum computing based on biomolecules . If $`1`$ and $`2`$ denote the two molecules, the coupling can be described by adding a term to the Hamiltonian

$`H_0`$ $`={\displaystyle \frac{J_{12}}{2}}\left(\sigma _+^1\sigma _{}^2+\sigma _{}^1\sigma _+^2\right)=J_{12}\left(\sigma _x^1\sigma _x^2+\sigma _y^1\sigma _y^2\right),`$ (20)
where $`\sigma ^\pm =\sigma _x\pm i\sigma _y`$. $`J_{12}`$ describes the FRET coupling and scales like $`1/R^3`$ where $`R`$ is the separation of the molecules .

What about the effect of the solvent? We assume that the molecules are sufficiently far apart that the reaction fields acting on the two biomolecules are uncorrelated. This assumption can be tested by looking at the wavevector dependence of the dielectric constant $`ϵ(\omega ,\stackrel{}{q})`$; this will be weak. We anticipate that the corresponding length scale is of the order of the separation of the polar molecules in the solvent which in turn is much larger than $`R`$, the separation of the biomolecules. Thus, we model the reaction fields by two independent collections of harmonic oscillators. The relevant Hamiltonian is then

$`H=`$ $`{\displaystyle \frac{1}{2}}ϵ_1\sigma _z^1+{\displaystyle \underset{\alpha }{}}\omega _\alpha a_\alpha ^{}a_\alpha +\sigma _z^1{\displaystyle \underset{\alpha }{}}g_\alpha (a_\alpha +a_\alpha ^{})+`$ (21)
$`{\displaystyle \frac{1}{2}}ϵ_2\sigma _z^2+{\displaystyle \underset{\beta }{}}\omega _\beta b_\beta ^{}b_\beta +\sigma _z^2{\displaystyle \underset{\beta }{}}\stackrel{~}{g}_\beta (b_\beta +b_\beta ^{})+`$
$`J_{12}(\sigma _x^1\sigma _x^2+\sigma _y^1\sigma _y^2),`$
and $`J(\omega ),\stackrel{~}{J}(\omega )`$ are as before, and $`ϵ_1,ϵ_2`$ are the energy gaps of the two biomolecules. In the Ohmic regime, this will have similarities to two Kondo impurities coupled by an XX interaction and in an external magnetic field . We have studied this problem in more detail in a separate article . There, we show that in most cases the above Hamiltonian may be reduced to a single spin-boson model. This allows us to give quantitative criteria, in terms of experimentally measurable system parameters, that are necessary for coherent Bloch oscillations of excitons between the chromophores, and suggest how these could be observed experimentally through FRET spectroscopy.

The above Hamiltonian (21) has a natural generalisation to many coupled biomolecules, such as the rings of chlorophyll molecules in a light harvesting complex in purple bacteria  or in a segment of base pairs in double stranded DNA . It will allow calculation of the localisation (the extent of spatial coherence) of excitons in such extended systems.

## 7 Summary

The major new result of this work is that we have shown that a chromophore and its environment may be modelled by an independent boson model, and derived the appropriate spectral density, which can be expressed in terms of experimentally obtainable parameters. We have estimated the critical coupling parameter $`\alpha `$, and found that quantum states are rapidly decohered by the interaction of the chromophore dipole with its surrounding polar environment. However, as can be seen from equation (13), the coupling to the environment is strongly dependent on the cavity radius, and further work must be done to determine the effects of surrounding proteins and solvation shells which might push the bulk solvent away from the chromophore and reduce the decoherence. Finally, we have introduced spin-boson models which could be used to model conformational change around a conical intersection and molecule coupled by fluorescence resonance energy transfer. These models, although simple, exhibit rich physics and may help to provide insight into the workings of these complex systems.

This work was supported by the Australian Research Council and by the University of Queensland Graduate School Research Travel Awards. We thank R.H. Austin, G.A.D. Briggs, I. Burghardt,P.C.W. Davies, B.L. Hu, B. Lovett, T. Martinez, P. Meredith,G.J. Milburn, J. Riesz and K. Schulten for helpful discussions. We thank the QIPIRC at Oxford for hospitality.