text:
# “Osiris”(HD209458b), an evaporating planet

## 1. Introduction

The extrasolar planet in the system HD209458 is the first one for which repeated transits across the stellar disk have been observed ($``$1.35% absorption; Henry et al., 2000; Charbonneau et al., 2000). Together with radial velocity measurements (Mazeh et al., 2000), this has led to a determination of the planet’s radius and mass, confirming that it must be a gas giant. The transits of this gaseous extrasolar planet offer a unique opportunity to investigate the spectral features of its atmosphere. Numerous searches for an atmospheric signature in the IR-optical wavelength range failed until the detection with Hubble Space Telescope (HST) of the dense lower atmosphere of HD209458b observed in the neutral sodium lines ($``$0.02% additional absorption; Charbonneau et al., 2002).

Far more abundant than any other species, hydrogen is well suited for searching atmospheric absorptions during transits, in particular over the stellar ultraviolet (UV) Lyman $`\alpha `$ emission line at 1215.67Å. With recent HST observations in the UV, Vidal–Madjar et al. (2003) detected the H i atomic hydrogen absorption ($``$15% ) over the stellar Lyman $`\alpha `$ line during three transits of HD209458b. Comparison with models showed that this absorption should take place beyond the Roche limit and could thus be understood only in terms of escaping hydrogen atoms.

The detection of this extended upper atmosphere (the exosphere) provides a new and crucial step in our understanding of the outer atmospheric structures of these “hot” extrasolar planets.

## 2. The observations

Three transits of HD209458b (noted A, B and C) were surveyed with the STIS spectrograph on-board HST ($``$20 km.s-1 resolution). For each transit, three consecutive HST orbits were scheduled such that the first orbit ended before the first contact to serve as a reference, the two following ones being partly or entirely within the transit. An average 15$`\pm `$4% (1$`\sigma `$) relative intensity drop (variation of the ratio between the flux in the “In” over the flux in the “Out” region, noted “In”/“Out”, see Fig. 1) near the center of the Lyman $`\alpha `$ line was observed during the transits. This is larger than expected for the atmosphere of a planet occulting only $``$1.5% of the star. Note that this drop has to be added on top of the $``$1.5% reduction of the whole stellar spectrum caused by the planetary disk transit which includes of course the Lyman $`\alpha `$ line.

Because of the small distance (8.5 R) between the planet and the star (allowing an intense heating of the planet and its classification as a “hot Jupiter”) the Roche lobe is only at 2.7 planetary radii (i.e. 3.6 RJup). Filling up this lobe with hydrogen atoms gives a maximum absorption of $``$10% during planetary transits. Since a more important absorption was detected, hydrogen atoms cover a larger area corresponding to a spherical object of 4.3 RJup. Observed beyond the Roche limit, these hydrogen atoms must be escaping the planet. Independently, the spectral absorption width (between $`130`$ and $`+100`$ km.s-1, see Fig. 1) also shows that some hydrogen atoms have large velocities relative to the planet, exceeding the escape velocity. This further confirms that hydrogen atoms must be escaping the planetary atmosphere.

### 2.1. Observational difficulties

These are essentially twofold : flux variations during transits could be either due to changes in the geocoronal Lyman $`\alpha `$ emission which covers and perturbs part of the stellar line or to intrinsic stellar flux variations.

First the geocoronal Lyman $`\alpha `$ emission produced by the hydrogen atoms present above HST in the Earth upper atmosphere (see e.g. Vidal–Madjar & Thomas, 1978) fills the aperture of the spectrograph. This results in an extended emission spread over the instrument detector, perpendicular to the dispersion direction. The extent of this emission along the slit allows however a precise correction at the position of the stellar spectrum. It was demonstrated that the geocoronal correction was made with high enough accuracy at least outside the central region ($`40`$ and $`+35`$ km.s-1, noted “Geo” in Fig. 1). No correlation was found between the geocoronal variations and the stellar Lyman $`\alpha `$ profile changes observed during transits showing independently that the corrections were properly made (see Fig. 2).

Second the observed Lyman $`\alpha `$ spectrum of HD209458 (G0V) is typical for a solar type star, with a double peaked emission originating in the stellar chromosphere. Stellar variability is thus a concern because it is well known that solar type stars present important flux variations in the far and extreme UV as it was observed in details in the case of the Sun (see e.g. Vidal–Madjar 1975). The total solar Lyman $`\alpha `$ flux varies by about a factor two, while its relative fluctuations in similar spectral domains (variations of the “In”/“Out” ratios) are smaller than $`\pm `$6%. Within a few months, a time comparable to the HD209458 observations, the “In”/“Out” solar flux ratio varies by less than $`\pm `$4% . If both stars behave in similar manners, this is an indication that the absorption detected is not of stellar origin. This is confirmed by the fact that no correlation was found between the total stellar flux changes (“In”+“Out” flux) and the absorption observed during transits (Fig. 2) showing that the observed variations are not related to stellar flux fluctuations.

## 3. The upper atmospheric model

The observed 15% intensity drop could only be explained if hydrogen atoms are able to reach the Roche lobe of the planet and then escape. To evaluate the amount of escaping atoms Vidal–Madjar et al. (2003) have built a particle simulation in which they assumed that hydrogen atoms are sensitive to the stellar gravity and radiation pressure (see Fig. 3). In their simulation, escaping hydrogen atoms expand in an asymmetric cometary like tail and are progressively ionized when moving away from the planet. Atoms in the evaporating coma and tail cover a large area of the star. A minimum escape flux of $`10^{10}`$ g.s-1 is compatible with the observations. However, due to saturation effects in the Lyman $`\alpha `$ absorption line, larger escape fluxes by several orders of magnitude would still be compatible with the observations.

To evaluate the actual escape flux, one has thus to estimate the vertical distribution of hydrogen atoms up to the Roche limit, in a more realistic heated up and tidally extended upper atmosphere.

### 3.1. The source of atomic hydrogen

Hi is produced through H2 photo–dissociation by extreme UV (EUV) radiation at the top of the atmosphere and through photolysis of hydrocarbons at the bottom. With an EUV heating source of 160 erg cm-2.s-1 on top of the atmosphere, in principle one can derive the amount of Hi atoms at the base of the upper atmosphere, in the region that could be considered as the source of atoms. Using realistic condition for the bottom of the atmosphere of HD209458b, Liang et al. (2003) have shown that the atomic hydrogen could be as high as 10% of H2 while in the Solar system giant planets it is of the order of 1% .This is one input parameter to be evaluated carefully.

### 3.2. The upper atmospheric temperature

Earlier estimates of the atmospheric escape flux concluded that the mass loss was not significant. However these estimates did not properly consider the temperature rise in the planetary upper atmosphere as it is observed in all Solar system planets. For example, the radiative equilibrium, effective temperatures Teff for the Earth and Jupiter are around 253 K and 150 K respectively while upper thermospheric and exospheric temperatures Tup are $``$1000 K for both planets (Chamberlain & Hunten 1987).

Burrows & Lunine (1995) estimated that the first “hot Jupiter” discovered, 51 Peg b, was comfortably within its Roche tidal lobe and thus that classical Jeans evaporation was insignificant. However they recognized that ultraviolet heating may play a role by mentioning that production of hot ions and atoms could be a more promising escape mechanism. By scaling from the case of Jupiter they evaluate that 1034 atoms.s-1 could possibly escape from 51 Peg b ($``$1010 g.s-1), a value similar to the Vidal–Madjar et al. (2003) lower limit.

Using the atmospheric effective temperature Teff=1300 K as the reference temperature, Guillot et al. (1996) concluded also that “hot Jupiters” even at such “high” temperatures are much too compact for classical Jeans escape. They also recognized that escape of hot ions is a possible and more likely process.

First Chassefière et al. (1996) considered upper atmospheric heating and evaluated that exospheric temperatures as high as Tup$``$104 K were possible leading to upper atmospheres extended up to 10 RJup. Such possibilities were again mentioned by Schneider et al. (1998) while Coustenis et al. (1998) underlined again that upper atmospheric temperatures of extrasolar planets should be quite different from the classical Teff evaluations. Mentioning upper atmospheric heating by stellar EUV, they evaluated Tup$``$4600 K in the case of 51 Peg b.

In the case of HD209458b, Seager & Sasselov (2000) made a detailed lower atmospheric model. They evaluated the lower atmosphere effective temperature to be Teff=1350 K. They did not consider any upper atmospheric heating. Following a parallel approach, Brown (2001) found Teff=1400 K. He however clearly stated that the model calculation did not involve He, H2 and H which can absorb UV photons and are very abundant, suggesting that an upper atmospheric heating could be expected.

Finally more recently Lecavelier des Etangs et al. (2004) evaluated the upper atmospheric temperature of HD209458b by taking into account FUV and EUV heating and different cooling mechanisms. They evaluated Tup$``$11000 K.

Note however that the observed high upper atmospheric temperatures in the Solar system are not all explained. In particular EUV solar fluxes are able to heat up the upper atmospheres of the giant planets by no more than few tens of degrees ($``$25 K in the case of Jupiter). Clearly other mechanisms as heating by solar wind particles and/or by gravity waves within the upper atmospheres may contribute (Hunten & Dessler 1977, McConnell et al. 1982). This suggests that the Lecavelier des Etangs et al. (2004) estimate could be a lower limit but knowing that more heating produces more escape and thus more cooling, more heating could have a small impact on the actual upper atmospheric temperature of the planet (see also Hébrard et al. 2004).

### 3.3. Turbulent transport

The transport of atomic hydrogen at high altitudes could be done via thermal diffusion but also through turbulent transport, as in the case of the Solar system giant planets. This transport is commonly described by the classical eddy diffusion coefficient parameter. It could vary by a huge 10000 factor, from $`10^4`$ cm2.s-1 at Uranus to $`10^6`$ cm2.s-1 at Jupiter to more than $`10^7`$ cm2.s-1 Saturn and finally $`10^8`$ cm2.s-1 at Neptune (see e.g. Atreya, 1986). Consequently, much more atomic hydrogen could be present at the source level (the lower thermosphere), and could considerably increase the amount of atomic hydrogen in the upper layers of the atmosphere. Eddy diffusion could thus change the vertical transport tremendously and possibly increase the amount of atomic hydrogen escaping the planet. In the case of the HD209458b, according to the probable extreme situation of that giant planet, one could expect such transport mechanisms. By ignoring the unconstrained turbulent transport, we will thus obtain a lower limit of the actual escape flux from the planet.

### 3.4. The atomic hydrogen escape flux

Estimating the source of atomic hydrogen, the upper atmospheric temperature, taking into account tidal forces and ignoring turbulent transports, Lecavelier des Etangs et al. (2004) were able to evaluate the vertical distribution of atoms in the upper atmosphere (Fig. 4). By extending the exobase concept to the level with no collisions up to the Roche limit (the exobase is usually defined by the level above which particles have no more collisions), they evaluated the number of atoms escaping the planet. They found that for upper atmosphere temperatures above $``$7000 K, the escape flux is larger than the flux derived from Lyman-$`\alpha `$ observations. But higher exospheric temperatures could be reached producing escape fluxes possibly as high as $``$1035 atoms.s-1 ($``$1011 g.s-1), i.e. about ten times larger than the minimum value evaluated from the observations. Such a high escape flux could put the escape into the blow-off regime in which all atmospheric species are carried away with atomic and molecular hydrogen (see e.g. Watson et al. 1981; Kasting & Pollack 1983; Chamberlain & Hunten 1987).

## 4. Lifetime of an evaporating planet

This raises the question of the lifetime of the evaporating extrasolar planets which may be comparable to stellar lifetimes. If so, the “hot Jupiters” could evolve faster than their star, eventually becoming “hot hydrogen–poor Neptune–mass planets”. The evaporation process, more efficient for planets close to their star might also explain the very few detections of “hot Jupiters” with orbiting periods shorter than 3 days (e.g., Konacki et al. 2003; Udry et al. 2003).

If the escape evaluation made by Lecavelier des Etangs et al. (2004) is correct, the planetary lifetimes could be evaluated as a function of both the star–planet distance and the planet initial mass (Fig. 4). This shows that very close to their stars, planets may have short lifetime. They could evolve into planets with hydrogen–poor atmospheres, or even with no atmosphere at all, their inner core being directly exposed.

## 5. Conclusion and naming planets

The unique case of the transiting planet HD209458b allowed the first exoplanets atmospheric studies. The case of HD209458b being so unique in the domain of exoplanets studies and its current name so uneasy to use (and unappealing to everyone outside the scientific community), it seems reasonable to propose another name in addition to the present ones already given by the standard nomenclature to its star (to be simply followed by the letter “b” to name the planet). The star names are as listed in the Centre de Données Stellaires (CDS): V\* V376 Peg, AG+18 2243, AGKR 19726, BD+18 4917, GSC 01688-01821, HD 209458, HIC 108859, HIP 108859, PPM 141002, SAO 107623, SKY# 42008, TYC 1688- 1821-1, uvby98 100209458, and YZ 18 9012.

The idea is to use the Egyptian God name “Osiris” who was killed and parts of his body spread over the whole Egypt. To have him come back to life his sister Isis searched for all the pieces. She found all of them but one and so, this God having lost part of his body, could perfectly represent HD209458b, the evaporating planet. This was proposed at the XIXth IAP Colloquium, Extrasolar Planets, Today & Tomorrow, on July 1, 2003.

## References

Atreya,K.A. 1986, in Atmospheres and Ionospheres of the Outer Planets and Their Satellites (Springer)

Brown,T.M. 2001, ApJ, 553, 1006

Burrows,A., & Lunine,J. 1995, Nature, 378, 333

Chamberlain,J.W., & Hunten,D.M. 1987, in Theory of planetary atmospheres (Orlando, Academic Press Inc.), International Geophysics Series 36, 1

Chassefière,E., et al. 1996, in Searching for Life in the Solar System and Beyond. A Research Discussion Meeting (The Royal Society, London)

Coustenis,A., et al. 1998, in ASP Conf. Ser., 134, 296

Charbonneau,D., Brown,T.M., Latham,D.W., Mayor,M. 2000, ApJ, 529, L45

Charbonneau,D., Brown,T.M., Noyes,R.W., Gilliland,R.L. 2002, ApJ, 568, 377

Désert,J.M., et al. 2004, these proceedings

Guillot,T., Burrows,A., Hubbard,W.B., et al., 1996, ApJ, 459, L35

Hébrard,G., et al. 2004, these proceedings

Henry,G.W., Marcy,G.W., Butler,R.P., & Vogt,S.S. 2000, ApJ, 529, L41

Hunten,D.M., & Dessler,A.J. 1977, Planet. Space Sci., 25, 817

Kasting,J.F., & Pollack,J.B. 1983, Icarus, 53, 479

Konacki,M., Torres,G., Jha,S., & Sasselov,D.D. 2003, Nature, 421, 507

Lecavelier des Etangs,A., Vidal–Madjar,A., et al., 2004, A&A, submitted

Liang,M.C., Parkinson,C.D., Lee,A.Y.T., & Yung,Y.L. 2003, astro-ph/0307037

Mazeh,T., et al. 2000, ApJ, 532, L55

McConnell,J.C., et al., 1982, Planet. Space Sci., 30, 151

Moutou,C., Coustenis,A., Schneider,J., et al., 2001, A&A, 371, 260

Sasselov,D.D. 2003, ApJ, 596, 1327

Schneider,J., et al. 1998, in ASP Conf. Ser., 134, 241

Seager,S., & Sasselov,D.D. 2000, ApJ, 537, 916

Udry,S., Mayor,M., Clausen,J.V., et al., 2003, A&A, 407, 679

Vidal–Madjar,A. 1975, Solar Physics, 40, 69

Vidal–Madjar,A., & Thomas,G.E. 1978, Planet. Space Sci., 26, 863

Vidal–Madjar,A., et al., 2003, Nature, 422, 143

Watson,A.J., Donahue,T.M., & Walker,J.C.G. 1981, Icarus, 48, 150


format: text

text:
# Solving the Richardson equations for Fermions

## Introduction

Exactly solvable models serve as a guideline for understanding the properties of correlated many-body systems. Already in the sixties, Richardson solved the eigenproblem of a constant pairing interaction in a set non-degenerate single-particle levels for Fermion Rich63  and Boson Richa  systems. However, it turned out that in the Fermion case the solutions exhibit singularities which are hard to treat numerically Rich64 ; Rich65 . Recently, exactly solvable pairing models have gained new attention Heri01 , with applications to nanometalic grains Sierra  (for a review see Delf01 ), Bose Einstein Condensates Duke01c  and nuclear physics Duke01b . It was first shown that the pairing model was integrable by finding the complete set of commuting integrals of motions Camb97 , and subsequently, three new families of fully integrable and exactly solvable models Duke01a  giving rise to a large class of pairing Hamiltonians with non-uniform matrix elements Dukec  were presented. These models are exactly solvable, except for the singularities occurring in Fermion systems for some critical values of the pairing strength. This problem, in spite of some early attempts to cure it Rich66a ; Ha , precluded for over forty years the use of these exactly solvable models for a wide range of applications, ranging from condensed matter to nuclear physics. Moreover, the recent developed extensions of the exact solution to pairing Hamiltonians including the isospin degree of freedom ($`O(5)`$ pairing) Pan  with promising applications to $`NZ`$ nuclei and high $`T_c`$ superconductivity, suffer from the same kind of singularities.

In this letter show how the Richardson equations can be solved numerically, avoiding the singularities, through an appropriate change of variables. This procedure provides a fast an accurate way to solve the equations for a constant pairing interaction. The method can be applied as well to more general exactly solvable Hamiltonians associated with a coupled set of non-linear equations of the Richardson type.

The exactly solvable pairing Hamiltonian has the following form:

$$H=\underset{j,m}{}e_ja_{jm}^{}a_{jm}\frac{g}{4}\underset{j,m,j^{},m^{}}{}a_{jm}^{}a_{j\overline{m}}^{}a_{j^{}\overline{m}^{}}a_{j^{}m^{}},$$
(1)
with $`e_j`$ any set of single-particle energies and $`g`$ the pairing interaction strength. The $`N`$-pair eigenstates of this Hamiltonian have the form

$$\underset{\alpha =1}{\overset{N}{}}\left[\underset{j,m}{}\frac{1}{2e_jx_\alpha }a_{jm}^{}a_{j\overline{m}}^{}\right]|0,$$
(2)
where $`|0`$ is a state without paired particles (a Racah quasispin vacuum state, Ring80 ). The corresponding energy $`E(\{x\})`$ is given by

$$E(\{x\})=0|H|0+\underset{\alpha }{}x_\alpha .$$
(3)
The complex variables $`x_\alpha `$ can be found by solving a set of non-linear equations:

$$\underset{j}{}\frac{d_j}{2e_jx_\alpha }+\underset{\beta =1,\beta \alpha }{\overset{N}{}}\frac{1}{x_\beta x_\alpha }+\frac{1}{2g}=0,$$
(4)
for $`\alpha =1,\mathrm{},N`$. The parameters $`d_j`$ depend on the level degeneracies and on the structure of the vacuum state $`|0`$. For Fermions, $`d_j=\frac{\nu _j\mathrm{\Omega }_j}{2},`$ where $`\mathrm{\Omega }_j`$ is the pair degeneracy (for the nuclear shell model $`\mathrm{\Omega }_j=j+1/2`$ ), and $`\nu _j`$ is the seniority of the level $`j`$. Note that $`d_j0`$ because of the Pauli principle. Solving this set of nonlinear equations solves the eigenproblem for the pairing Hamiltonian Eq.(1). Unfortunately, the algebraic solution becomes numerically unstable at certain critical values of the interaction strength Rich64 ; Rich65 . This is caused by singularities in the first and second terms in Eq.(4), when some of the variables $`x_\alpha `$ are approaching the value $`2e_j`$. One can understand this from the electromagnetic analogy for the exactly solvable pairing model Duke02 : in the Fermion case, the single-particle levels $`e_j`$ and the variables $`x_\alpha `$ correspond to opposite charges. Therefore a group of variables can cluster around a single-particle level in such a way that for each of the variables the repulsive charge of the other variables is compensated for by the attractive charge of the single-particle level. These singularities occur for Fermions in double or multiple degenerate levels. In the case of doubly degenerate equidistant levels these singularities can be handled for the ground state Rich66a , but the problem arises again in the treatment of the excited states Roman . In previous calculations it was necessary to tune the solutions by hand as soon as they approach a singularity. No general solution method was known until now. Figure 1 illustrates the behavior of the variables in case of multiply degenerate levels (see below for the details of the model).

A general approach to solve the Richardson equations Eq.(4), starts from an approximate solution in the weak-interaction limit (see below). Then this solution is evolved adiabatically up to the desired interaction strength by gradually increasing the value of $`g`$. At each step in $`g`$, the previous solution has to be updated. For dealing with the singularities, it is useful to notice first that the variables $`x_\alpha `$ in Eq.(4) are most sensitive to the other variables nearby. Therefore one can divide the set of variables $`x_\alpha `$ into several clusters of variables, grouped around different single-particle levels. By solving the equations for each cluster separately, one can obtain a solution for the whole system iteratively. A practical way to organize the clusters, is to link each variable to its nearest single-particle level $`2e_j`$, and to consider a cluster for each level that has variables around it.

The question now is how to solve the equations for each cluster, particularly in the case of singularities. Let us consider the set of indices $`C_k`$ of the $`N_k`$ variables that cluster around a level $`2e_k`$. Then one can consider the equations

$$\frac{d_k}{2e_kx_\alpha }+\underset{\beta C_k,\beta \alpha }{}\frac{1}{x_\beta x_\alpha }+F_k(x_\alpha )=0,\alpha C_k,$$
(5)
with

$$F_k(x)=\frac{1}{2g}+\underset{j,jk}{}\frac{d_j}{2e_jx}+\underset{\beta ,\beta C_k}{}\frac{1}{x_\beta x}.$$
(6)
The function $`F_k(x)`$ describes the influence of the other levels and the variables of the other clusters on the variables in the cluster $`C_k`$. Because of the way the clusters are set up, the function $`F_k(x)`$ will be a smooth function in the region around $`e_k`$ where the variables of the cluster are located. The singularities will occur in the first two terms of Eq.(5). In the case that some of the variables in the cluster approach the value $`2e_k`$, the divergences in the first and the second term of Eq.(5) must cancel out. Multiplying Eq.(5) by $`2e_kx_\alpha `$, and summing over the $`n_k`$ variables $`x_\alpha `$ at the singular point, leads to the condition

$$n_k=2d_k+1,$$
(7)
with $`n_k`$ the number of variables that actually converge to $`2e_k`$. For Fermions the value $`(2d_k)`$ corresponds to the pair degeneracy of the level $`e_k`$: because of the Pauli principle, no more Fermion pairs can occupy the level. Trying to put more pairs in that level results in a singularity. In fact, the structure of the ground state Eq.(2) does not result in a forbidden occupation of the level. However, on expanding the wave function of Eq.(2) in terms of the pair creation operators $`a_{jm}^{}a_{j\overline{m}}^{}`$, one finds that the leading term cancels out because of the Fermionic anti-commutation rules. This translates into the numerical difficulties encountered in the solution of the equations. For Bosons $`d_k>0`$, and hence Eq.(7) shows that singularities do not occur in the Bosonic case. This can also be understood from the electrostatic analogy: Boson pairs and single-particle levels have charges of the same sign. Therefore the variables try to avoid each other at all times, and singularities do not occur.

The above procedure suggests a way to remove the singularities from the equations: multiplying Eq.(5) with $`(2e_kx_\alpha )^p`$, for some power $`p`$, and summing over all variables $`x_\alpha `$ in the cluster. The resulting equations become, for $`p>1`$:

$$\left(d_k+N_k\frac{p}{2}\right)S_{p1}+\frac{1}{2}\underset{k=2}{\overset{p1}{}}S_{k1}S_{pk}+R_p=0,$$
(8)
with

$`S_p`$ $`=`$ $`{\displaystyle \underset{\alpha C_k}{}}(2e_kx_\alpha )^p`$ (9)
$`R_p`$ $`=`$ $`{\displaystyle \underset{\alpha C_k}{}}(2e_kx_\alpha )^pF_k(x_\alpha ).`$ (10)
The compact form of Eq.(8) suggests that it might be advantageous to solve them for the new variables $`S_p`$ instead of the original variables $`x_\alpha `$. Note that given a set of variables $`S_1,\mathrm{},S_{N_k}`$, one can easily construct the polynomial whose roots correspond to the values $`2e_kx_\alpha `$. Hence one can switch from one set of variables to the other. The problem comes with the quantities $`R_p`$: these are functions of the $`x_\alpha `$, and it is not straightforward to express them as functions of the $`S_p`$. However, for a given set of variables $`x_\alpha `$, one can easily evaluate the values $`S_p`$ and $`R_p`$. Furthermore, one can evaluate the gradient matrix $`G`$, with $`G_{lm}`$ the derivative with respect to $`S_m`$ of Eq.(8) for $`p=l+1`$,

$$G=G^S+G^R,$$
(11)
where

$`G_{lm}^S`$ $`=`$ $`\{\begin{array}{cc}d_k+N_k\frac{l+1}{2},\hfill & \text{for }m=l\hfill \\ S_{lm},\hfill & \text{for }m<l\hfill \\ 0,\hfill & \text{for }m>l\hfill \end{array}`$ (15)
$`G_{lm}^R`$ $`=`$ $`{\displaystyle \frac{dR_{l+1}}{dS_m}}={\displaystyle \underset{\alpha C_k}{}}{\displaystyle \frac{dR_{l+1}}{dx_\alpha }}{\displaystyle \frac{dx_\alpha }{dS_m}},`$ (16)
for $`l,m=1,\mathrm{},N_k`$. $`G^R`$ can be evaluated accurately using a special inversion algorithm for Vandermonde matrices Golu89 . Therefore one can solve the new set of equations, Eq.(8) for $`p=2,\mathrm{}N_k+1`$, in the new variables $`S_1,\mathrm{}S_{N_k}`$ using a standard gradient technique such as the multi-dimensional Newton-Raphson method Pres92 . However, in the case of a singularity, the gradient matrix becomes ill-conditioned: the diagonal elements of the gradient matrix are given by

$$G_{ll}^S=d_k+N_k\frac{l+1}{2},$$
(17)
for $`l=1,\mathrm{}N_k`$. The diagonal element will vanish for the index $`l_s=2(d_k+N_k)1`$. This will occur as soon as $`N_k2d_k+1`$, which matches the value for which singularities can occur, see Eq.(7). In such a case the lower-triangular matrix $`G^S`$ becomes singular. The other part of the gradient matrix, $`G^R`$, is derived from the smooth function $`F_k`$. A series expansion of $`F_k(x)`$ in $`x`$ will be dominated by the lowest orders. Therefore the elements of $`G_{lm}^R`$ are very small for larger values of $`m`$. As a result, the value of $`S_{l_s}`$ can not be determined accurately from the set of equations Eq.(8). One can avoid this problem by limiting the cluster sizes to at most the critical value $`N_k=n_k=2d_k+1`$, and by using $`gS_1`$ as an unknown variable instead of $`S_{N_k}`$. If more variables are found near to the same single-particle level, one can always divide the cluster into smaller, well-separated clusters, because at most $`n_k`$ of the variables can approach the single-particle level closely. Knowing $`S_1`$ and $`S_1,\mathrm{},S_{N_k1}`$, one can still straightforwardly construct the polynomial whose roots give the corresponding values $`x_\alpha `$. Therefore one can easily switch between the two sets of variables. Furthermore $`gS_1`$ behaves smoothly, even at a singularity. To set up an efficient gradient method, it is useful to replace Eq.(8) for the last value, $`p=N_k+1`$, by a similar equation obtained using $`p=0`$:

$$d_kgS_1+gR_0=0,$$
(18)
with $`R_0=_{\alpha C_k}F_k(x_\alpha )`$.

For weak interaction strengths the function $`F_k(x_\alpha )`$ is dominated by the constant term $`\frac{1}{2g}`$, see Eq.(6). In the weak interaction limit one can take $`F_k`$ to be a constant. The resulting functions $`R_p`$ take the simple form

$$R_p=\frac{S_p}{2g}.$$
(19)
Now the equations Eq.(8) and Eq.(18) can be solved straightforwardly to yield the variables $`S_l`$, from which one can construct the polynomial that gives a unique set of variables $`x_\alpha `$. The resulting eigenstate will depend on the size of the cluster for each of the single-particle levels. This establishes a one-to-one correlation between the eigenstates of the non-interacting system ($`g=0`$) and the eigenstates of the weakly-interacting system. One can conclude that the Richardson equations are complete: their solutions generate all eigenstates, and there are no spurious solutions.

One can obtain the solutions for strong interaction strengths by solving the weakly interacting case first, and then gradually increasing the interaction strength. At each step one can use the gradient method outlined above in order to update the solution to the new interaction strength. It is useful to adapt the stepsize in interaction strength to the convergence of the iterative procedure by taking smaller steps in $`g`$ when the convergence of the Newton-Raphson method for the variables $`S_i`$ becomes slower, which typically occurs around the critical $`g`$ values. One more ingredient is needed to avoid problems with the singularities: when the interaction strength passes through a critical value, the variables $`x_\alpha `$ passing through a singularity can change from real to complex or vice versa. At the same time the variables $`S_l`$ will become very small, except for $`l=1`$. Even using the new variables, the gradient method does not lead to the right solution when it has to pass through critical values of the interaction strength. One can avoid this problem by including an extrapolation step based upon the previous solutions. This extrapolation has to be done in the variables $`S_l`$, because they vary smoothly through the singularities and they remain real all the time. Assume that converged solutions $`x_\alpha ^{}`$ and $`x_\alpha ^{\prime \prime }`$ were obtained for values $`g^{}`$ and $`g^{\prime \prime }`$ of the interaction strength. The variables $`x_\alpha ^{\prime \prime }`$ are grouped into clusters. For each cluster, one evaluates the variables $`S_p^{}`$ and $`S_p^{\prime \prime }`$. Because the variables $`S_p`$ behave smoothly, even near a singularity, one can estimate the variables $`S_p`$ for the new interaction strength $`g`$ by linear extrapolation:

$$S_p=\frac{(gg^{})S_p^{\prime \prime }(gg^{\prime \prime })S_p^{}}{g^{\prime \prime }g^{}}$$
(20)
The resulting values of $`S_p`$ can then be updated using the Newton-Raphson method or another gradient method. The extrapolation step avoids the problems with the singularities and greatly improves the convergence of the method.

As an example, consider the level scheme listed in Table 1, together with a constant pairing interaction. This model describes neutrons in 56Fe; its ground state and finite-temperature properties have been studied using a Quantum Monte Carlo method Romb98 . The eigenstates can also be found through the solution of Richardson’s equations, Eq.(4), or through Lanczos diagonalization in a seniority basis Voly01 . The full many-body space has a dimension of the order of $`10^{15}`$, while the zero-seniority basis has dimension 14894. In Fig. 2 the lowest zero-seniority eigenvalues for this model are shown as a function of the interaction strength, calculated using the Lanczos method and using the method explained above (only 50 Lanczos iterations were used).

To calculate the ground state, a straightforward implementation of the Newton-Raphson method for the original equations, Eq.(4), works well up to an interaction strength of $`g0.2`$. It turns out that a singularity occurs around the $`1d_{3/2}`$ level at a value of $`g=0.245`$. A Fortran-95 computer program was written based upon the procedure outlined above. It was able to solve the equations for all interaction strengths in a matter of seconds. Figure 3 shows the behavior of the three $`x`$ variables that cluster around the $`1d_{3/2}`$ level as a function of $`g`$. In Fig. 4 one can see that the corresponding variables $`S_l`$ behave much more smoothly. More singularities occur around other levels at higher interaction strengths, as is shown in Fig. 1. The solution of the Richardson equations is faster than the Lanczos method, it requires less computer memory, and it works for all the eigenstates, not just the ground state. Moreover, the procedure can deal with much larger systems, well beyond the limits of large-scale exact diagonalizations.

This work shows that the exactly solvable pairing models are indeed solvable in practice, even for Fermions with multiple degeneracies. It opens up a whole new range of applications for these models. Furthermore, similar techniques might be useful to solve the non-linear equations for other exactly solvable pairing and spin models Duke01a ; Feng98 ; Pan ; Shas

###### Acknowledgements.

We wish to thank R. Richardson, S. Pittel, J. Draayer, F. Pan and K. Heyde for the interesting discussions and suggestions. This work was supported by the Fund for Scientific Research - Flanders (Belgium), the Research Board of Ghent University and by the Spanish DGI Grant BFM2003-05316-C02-02.


text: # Untitled Document


text:
# Loop Quantum Cosmology and Boundary Proposals

## I Introduction

Loop quantum cosmology is the cosmological sector of loop quantum gravity,Rov:Loops ; ThomasRev  a canonical quantization of general relativity using Ashtekar’s variables. In this formulation, gravity becomes a constrained gauge theory with the basic classical fields being an SU(2)-connection $`A_a^i`$, a sum of the spin connection and the extrinsic curvature, and its conjugate, the densitized triad $`E_i^a`$. These fields are represented in the quantum theory by holonomies of the connection and fluxes of the triad, which allows a mathematically well-defines formulation of the full theory.ALMMT  Nevertheless, the complexity of the full formalism makes it advantageous to consider simplified models in a minisuperspace quantization. Since there is a well-developed mathematical structure behind the full theory, the reduction to symmetric models can be done in such a way that essential features of the full theory, e.g., the discreteness of geometric spectra, are preserved. This is in contrast to the older Wheeler–DeWitt quantizationDeWitt ; Misner ; QCReview  where the full theory would not even be known beyond a purely formal level.

Technically, states of the full theory are functions on the space of connections which depend on $`A_a^i`$ only via holonomies (cylindrical functions). Symmetric states are by definition distributions on the space of connections which are supported only on the sub-manifold of invariant connections.SymmRed  It turns out that such a reduction even to the simplest cosmological models yields a quantum representation which is inequivalent to the one chosen in the Wheeler–DeWitt quantization (which is possible also in systems with a finite number of degrees of freedom since the Weyl algebra is not represented continuously).Bohr  This translates into a different point of view as to the representation of basic variables. While the Wheeler–DeWitt quantization assumes that the extrinsic curvature can be promoted to a well-defined operator, the background independent techniques employed in the loop quantization allow only a representation of holonomies as operators, not the connection itself which in the isotropic setting corresponds to the extrinsic curvature. Also triad operators behave differently: they would have continuous spectra in a Wheeler–DeWitt quantization while they have discrete spectra in a loop quantization. On the other hand, these basic properties are the same in full loop quantum gravity and in loop quantum cosmology such that one can regard those cosmological models as reliable tests.

These basic technical details have far-reaching physical consequences. For instance, the Wheeler–DeWitt quantization cannot eliminate the cosmological singularities that plague general relativity, while they do not occur in loop quantum cosmology.Sing  Also the issue of boundary conditions, which has been the main issue in the Wheeler–DeWitt quantization, appears in a different light.DynIn ; Essay  All these effects are intimately linked to the underlying discreteness of geometry.

In this contribution we describe both the basic formulation of the theory, emphasizing differences to the Wheeler–DeWitt quantization, and physical applications. Both parts can be read independently of each other.

## II Basic Formalism of Loop Quantum Gravity and Loop Quantum Cosmology

Loop quantum gravity is based on general relativity in Ashtekar’s canonical variablesAshVar ; AshVarReell  $`A_a^i=\mathrm{\Gamma }_a^i(E)\gamma K_a^i`$ and $`E_i^a`$ on a spacelike slice $`\mathrm{\Sigma }`$, where $`\mathrm{\Gamma }_a^i`$ is the spin connection associated with the densitized triad $`E_i^a`$, $`K_a^i`$ the extrinsic curvature, and $`\gamma `$ the Barbero–Immirzi parameter which is a positive real number. For a background independent quantization it is important to represent these basic fields via holonomies

$$h_e(A)=𝒫\mathrm{exp}_eA_a^i\tau _i\dot{e}^a𝑑t$$

where $`e`$ is a curve in $`\mathrm{\Sigma }`$ with tangent vector field $`\dot{e}^a`$ and $`\tau _i`$ are Pauli matrices, and fluxes

$$F_S^{(f)}(E)=_SE_i^af^in_ad^2y$$

where $`S`$ is a surface in $`\mathrm{\Sigma }`$ with co-normal $`n_a`$ and $`f^i`$ a test function. If all curves, surfaces and test functions are allowed, these quantities have the full information contained in the fields $`A_a^i`$ and $`E_i^a`$.

The variables are constrained by the usual SU(2)-Gauss constraint as well as the diffeomorphism constraint, which requires invariance under deformations of space, and the Hamiltonian constraint which yields the dynamics of the theory. Since the diffeomorphisms act directly on the basic fields, a representation of their classical Poisson-\*-algebra has to be diffeomorphism invariant if this classical symmetry is to be preserved in the quantization.

### II.1 Representation of the Basic Classical Algebra

There is in fact a natural diffeomorphism invariant representation of the holonomy/flux algebra which has recently been shown to be unique.FluxAlg ; Meas ; SuperSel ; Irred ; HolFluxRep  A heuristic construction of this representation can be done as follows. Since there is a rich mathematical structure on the space of connections, it is convenient to look for the state space in the connection representation. The simplest such state is the “ground state” 1$`(A)=1`$ which does not depend on the connection at all. All higher states can then be “created” by acting with holonomies as multiplication operators. Since one can act with several holonomies, one at a time, the final states will depend on the connection non-trivially, but only along the edges in the holonomies chosen, and only via the holonomies (cylindrical functions). Labels of a given state, the edges of the holonomies used to construct it, can be arranged to a graph in $`\mathrm{\Sigma }`$ which may have self-intersections and also knotting and linking information. Furthermore, since one can act with the same holonomy several times, one also needs a label representing the number of times a given edge occurs in the graph. Another way to include this information is by using irreducible SU(2)-representations since acting once with a holonomy creates dependence of the state via the fundamental representation while acting several times then creates dependence according to the SU(2)-coupling rules. An independent set of states can be seen to be given by the spin network basisRS:Spinnet  given by graphs whose edges are labeled by irreducible SU(2)-representations, or half-integer spin labels $`j`$.

While any such state depends on the connection only via finitely many holonomies, the full state space is obtained as a projective limit over the set of graphs.FuncInt ; ALMMT  Thus, there is no restriction on the dependence of states on connections and one obtains a faithful representation of the infinite dimensional space of connections. This state space carries a natural diffeomorphism invariant inner product, the Ashtekar–Lewandowski inner product, which defines the Hilbert space of loop quantum gravity.

Fluxes represent the conjugate of the connection and thus act as (functional) derivative operators in the connection representation. Since the states depend on the connection via holonomies, the action of flux operators can be obtained from the chain rule. A flux operator associated with a surface $`S`$ has a non-trivial action on a state labeled by a graph $`g`$ only if $`S`$ and $`g`$ intersect each other transversally. If they do, the action is given by a sum over transversal intersection points, each contribution counting the multiplicity of the edge encoded in its spin label $`j`$. Thus, flux operators count the intersection number of the surface with the graph of a state, counting also multiplicities, which is an integer. In this way, flux operators obtain a discrete spectrum which translates to discrete spectra of other geometrical operators constructed from the densitized triad, such as area or volume.AreaVol ; Area ; Vol2

### II.2 Cosmological Models

In order to follow the reduction procedure, we need to know the sub-manifold of invariant connections for an isotropic symmetry group. Any isotropic connection can be written in the formcosmoI

$$A_a^i=c\mathrm{\Lambda }_I^i\omega _a^I$$

where $`\omega ^I`$ are left-invariant 1-forms on $`\mathrm{\Sigma }`$ for the action of the symmetry group, $`\mathrm{\Lambda }`$ is an internal triad (an $`SO(3)`$-matrix) which is purely gauge, and $`c`$ is the only gauge-invariant component of the connection. Similarly, an isotropic triad has the dual form

$$E_i^a=p\mathrm{\Lambda }_i^IX_I^a$$

with left-invariant vector fields $`X_I`$ dual to $`\omega ^I`$.

Isotropic states are then distributions in the full state space which are supported only on isotropic connections. Thus, gauge invariant isotropic states can depend only on the component $`c`$, and an explicit reduction, e.g. by restricting spin network states to isotropic connections, shows that any isotropic state can be represented as

$$|\psi =\underset{\mu }{}\psi _\mu |\mu $$

where the sum is over a countable subset of $``$. The isotropic states

$$c|\mu =e^{i\mu c}$$
(1)
form an orthonormal basis since their inner product is given by

$$\mu |\mu ^{}=\delta _{\mu ,\mu ^{}}$$
(2)
which implies that the Hilbert space is non-separable (all $`\mu `$ are allowed).

In this connection representation, holonomies $`e^{i\mu c}`$ act as multiplication operators

$$e^{i\mu c}|\mu ^{}=|\mu +\mu ^{}$$
(3)
and the flux $`p`$ as a derivative operator

$$\widehat{p}|\mu =\frac{1}{6}\gamma \mathrm{}_\mathrm{P}^2\mu |\mu .$$
(4)
While the Hilbert space and the action of basic operators are much simpler than in the full theory, their technical features are the same: As already noted, the Hilbert space is non-separable. Furthermore, only holonomies, not $`c`$ itself are promoted to well-defined operators; it is impossible to derive an operator for $`c`$ from the holonomy operators by differentiation since the operator family $`(e^{itc})_t`$ is not continuous in $`t`$ at $`t=0`$:

$$\mu |e^{itc}|\mu =\mu |\mu +t=\delta _{t,0}.$$
(5)
Also as in the full theory, the flux operator $`\widehat{p}`$ has a discrete spectrum: its eigenstates $`|\mu `$ are normalizable (even though a continuous range of eigenvalues is allowed).

The results (1), (2), (3) and (4) can also be used as axioms to define isotropic loop quantum cosmology, without following the reduction from the full theory. This may then be seen as a quantization alternative to and, as we will see, different from the Wheeler–DeWitt quantization (see also \[BohrADM, \] for a similar discussion in terms of ADM variables).

In order to compare this basic framework with that of the Wheeler–DeWitt quantization, we start with the formal Hilbert spaces. In a Wheeler–DeWitt quantization one would simply employ a Schrödinger quantization of the single degree of freedom $`c`$ such that the Hilbert space is $`L^2(,dc)`$ and the conjugate momentum $`p`$ acts as a derivative operator. The loop quantization is based on a different, inequivalent Hilbert space since it is non-separable. It can be written as $`L^2(\overline{}_{\mathrm{Bohr}},d\lambda )`$ where $`\overline{}_{\mathrm{Bohr}}`$ is the Bohr compactification of the real line $``$. It is a compact Abelian group which contains $`\overline{}_{\mathrm{Bohr}}`$ as a dense subset. The measure $`d\lambda `$ of the isotropic loop Hilbert space is the Haar measure on this compact group.

Thus, the Wheeler–DeWitt quantization and the loop quantization are based on inequivalent Hilbert spaces such that they cannot be isomorphic to each other. This is possible even with finitely many degrees of freedom since the Weyl algebra $`(e^{itc},e^{isp})_{t,s}`$ is not represented continuously in the loop representation; see (5). This assumption of the Stone-von Neumann theorem, which is usually used to prove the uniqueness of the quantum mechanical representation, is thus violated. As a consequence we have the basic differences for the operators: in the Wheeler–DeWitt quantization there is an operator $`\widehat{c}`$, which does not exist in the loop case, and the operator $`\widehat{p}`$ of the Wheeler–DeWitt quantization has a continuous spectrum compared to the discrete one in the loop quantization.

On the other hand, the loop quantization of the isotropic model shares basic features with full loop quantum gravity. In the full theory, the quantum configuration space is a compactification $`\overline{𝒜}`$ of the classical space of connections $`𝒜`$ which is contained in $`\overline{𝒜}`$ as a dense subset. This can be seen by analyzing the projective limit construction sketched above.ALMMT  The situation in the isotropic case is completely analogous: the isotropic space of connections $`𝒜_{\mathrm{iso}}`$ is extended and compactified to the quantum configuration space $`\overline{𝒜}_{\mathrm{iso}}\overline{}_{\mathrm{Bohr}}`$ which contains $`𝒜_{\mathrm{iso}}`$ as a dense subset. Here, the same mechanism occurs just in a simpler context. As consequences we have both in the full theory and in the isotropic case a non-separable Hilbert space, no well-defined operators for connections, and flux operators with discrete spectra. There is a difference, though, namely that the isotropic flux operator has a continuous set of eigenvalues even though its spectrum is discrete. This property comes from a degeneracy manifest in the reduction procedure: the label $`\mu `$ of an isotropic state $`e^{i\mu c}`$ does not distinguish between the parameter length of an edge and the spin when we compare it to a holonomy. There are technical subtleties resulting from this issue that can be dealt with appropriately (see \[Bohr, \] for a detailed discussion of all these concepts). Isotropic models, then, provide reliable tests of basic and important features of the full theory. Calculations in these models are certainly much simpler such that explicit results can be obtained. This is already shown by the volume operator which is rather complicated in the full theory while its spectrum can be obtained explicitly in the isotropic context.cosmoII  Since it plays an important role in constructing other operators such as the Hamiltonian constraint,QSDI ; QSDV  also the dynamical equations become more treatable. We now turn to the physical applications which can be derived in this way.

## III Cosmology

The dynamical equation of isotropic cosmology is the Friedmann equation which is the same as the Hamiltonian constraint equation. In connection variables, the isotropic gravitational part of the constraint reads

$$H=12\gamma ^2\kappa ^1(c(ck)+(1+\gamma ^2)k^2/4)\sqrt{|p|}$$
(6)
with the gravitational constant $`\kappa =8\pi G`$ and the curvature parameter $`k`$ which is zero for the spatially flat model and one for the closed model. In order to relate this to the usual form of the Friedmann equation, we need to know the relation between the basic variables $`(c,p)`$ and the scale factor $`a`$. For the triad component, one has $`|p|=a^2`$ ($`p`$ can be negative, indicating the orientation of the triad), while the connection component is $`c=(k\gamma \dot{a})/2`$. If we insert this into the constraint equation $`H+H_{\mathrm{matter}}=0`$ with some matter Hamiltonian, we obtain in fact the Friedmann equation

$$3(\dot{a}^2+k^2)a=8\pi Ga^3\rho (a).$$

### III.1 Difference Equation

The quantization of the Hamiltonian constraint is rather complicated (see \[Closed, \] for the expression) since it is a composite operator in terms of the basic ones, holonomies and fluxes. To solve the constraint, we need to act on a general state $`|\psi =_\mu \psi _\mu (\varphi )|\mu `$ where the coefficients depend on matter fields which we collectively denote as $`\varphi `$. It is more convenient to write down the constraint equation for the coefficients $`\psi _\mu (\varphi )`$ directly, i.e. to work in a triad rather than connection representation, which we can later compare to the Wheeler–DeWitt constraint. In this way, our dynamical law becomes the difference equation

$`(V_{\mu +5}V_{\mu +3})e^{ik}\psi _{\mu +4}(\varphi )(2+\gamma ^2k^2)(V_{\mu +1}V_{\mu 1})\psi _\mu (\varphi )`$
$`+(V_{\mu 3}V_{\mu 5})e^{ik}\psi _{\mu 4}(\varphi )={\displaystyle \frac{1}{6}}\gamma ^3\kappa \mathrm{}_\mathrm{P}^2\widehat{H}_{\mathrm{matter}}(\mu )\psi _\mu (\varphi )`$ (7)
using the volume eigenvalues

$$V_\mu =(\gamma \mathrm{}_\mathrm{P}^2|\mu |/6)^{3/2}.$$
(8)
A Wheeler–DeWitt quantization, in contrast, would result in a differential equation which turns out to be a good approximation to (III.1) at large volume, $`|\mu |1`$, where the differences can be Taylor expanded.SemiClass  At small $`\mu `$, however, the discrete formulation is very different from the continuous one. This is crucial when we consider the fate of the classical singularity which is not removed in the Wheeler–DeWitt quantization: generically, curvatures and energy densities still diverge and the classical singularity presents a boundary to the evolution (in some special models one can argue for a more regular behavior in a different sense, see e.g. \[Coule, \]). In the loop quantization, on the other hand, it turns out that there is no singularity;Sing ; IsoCosmo  instead, the dynamical law provides a unique evolution from the wave function at positive $`\mu `$ to its values at negative $`\mu `$. This can be seen directly when one starts with initial data at large positive $`\mu `$ and uses (III.1) as a recurrence relation to evolve backward toward the classical singularity at $`\mu =0`$. There are no problems at first as long as the lowest order coefficient, $`V_{\mu 3}V_{\mu 5}`$, in the difference equation is non-zero. For $`\mu =4`$, however, this coefficient is zero which means that $`\psi _0(\varphi )`$ remains undetermined in this way. If this value would then appear in the next stages of the recurrence, the evolution would break down and the classical singularity would present a boundary also for the quantum evolution. As one can see, this does not occur: also for $`\mu =0`$ and $`\mu =4`$, where $`\psi _0`$ would appear in (III.1) at first sight, it drops out thanks to vanishing coefficients. Here we use the fact that quantum geometry Hamiltonians imply the special featureInvScale  $`\widehat{H}_{\mathrm{matter}}(0)=0`$ which is independent of quantization ambiguitiesAmbig  and comes from a general procedure in the full theory.QSDV  Thus, $`\psi _0`$ also drops out in the matter part of the constraint equation and we can follow the evolution completely to negative $`\mu `$. (For a classical picture of the corresponding bounce in a closed model, see \[BounceClosed, \].)

The classical singularity, therefore, does not appear as a natural starting point where one could pose initial conditions. Nevertheless, the issue of boundary conditions appears in a new light now.DynIn ; Essay  When following the evolution we noticed that $`\psi _0`$ remains undetermined and we ignored the part of the constraint where it dropped out the first time ($`\mu =4`$). However, this part of the constraint still has to be satisfied providing us with a linear relation between $`\psi _4`$ and $`\psi _8`$. This, in turn, translates to a linear relation between the initial data we chose. In the isotropic case this is just the right amount to specify the gravitational wave function completely up to its norm (ignoring the matter field dependence) as usually done in the Wheeler-DeWitt quantization. In order to compare with Wheeler–DeWitt wave functions we have to assume that the solution we obtain is almost constant at large volume at least for positive $`\mu `$, i.e. that it does not have wild oscillations at the Planck scale. This condition (pre-classicalityDynIn ; FundamentalDisc ) ensures that the solution has a continuum limit to be compared with a Wheeler–DeWitt wave function. There are two independent such solutions which, when using the linear relation, are reduced to a unique (up to norm) linear combination. In this way, the dynamical law implies the boundary condition (dynamical initial condition) which is not chosen separately.

### III.2 Boundary Proposals

Boundary proposals for the Wheeler–DeWitt quantization have been studied for some time, mainly focusing on the tunnelingtunneling  and no-boundary proposals.nobound  We now have a very different formulation in loop quantum cosmology with a discrete evolution equation and derived boundary conditions. Nevertheless, after evolving the wave function to large volume, where the Wheeler–DeWitt quantization is a very good approximation, one can compare different boundary proposals and possible physical implications. For this purpose, we will only consider positive $`\mu `$ and $`p`$ for the rest of this section.

The first proposal is due to DeWitt who suggested that the wave function should vanish at the classical singularity,DeWitt  $`\psi (0)=0`$. The dynamical initial conditions in fact coincide with this one for simple systems such as de Sitter space with flat slices (see Fig. 1).DynIn  In more complicated cases, however, the discreteness at small volume implies differences in the behavior of the wave functions. In this way, the dynamical initial conditions even work for models where DeWitt’s condition would not be well-posed.Scalar

Other boundary proposals have almost exclusively been studied for the closed model, which is in particular true for the tunneling and no-boundary proposals. In the formulation of the proposals, one simply considers the closed model with a cosmological constant, where there is a classically forbidden region between values of the scale factor between zero and some minimum classically allowed value. There are detailed evaluations of their different implications, but the central difference we are concerned with here is as follows. The tunneling wave function starts with a large value at $`a=0`$ and then decays exponentially in the classically forbidden region, which is motivated as a tunneling process of the universe. After that, the wave function becomes oscillating in the classically allowed region. The no-boundary proposal, on the other hand, does have a non-vanishing exponentially increasing contribution to the wave function which will soon dominate. Thus, the wave function starts with a small value at $`a=0`$ and grows much higher before reaching the oscillatory phase. This different behavior is claimed to have consequences for chaotic inflation when one computes the probability distribution for the initial inflaton. This issue is still disputed, but it is often said that the no-boundary proposal prefers small initial values for the inflaton which would result in an insufficient amount of inflation.

In order to see how the loop wave function fits into this situation it is sufficient to note that the tunneling wave function is very special. It requires that the exponentially growing part does not contribute to the wave function which can be achieved only in this one very special case. Any other proposal would have a non-vanishing contribution of the growing mode and would soon differ from the tunneling wave function. Thus, any proposal which is not tailored to be exactly the tunneling one, would be much closer to the no-boundary condition at large volume. This is in particular true for the loop wave function with its dynamical initial condition. A numerical analysis (Figs. 2 and 3) shows that this wave function also increases exponentially in the classically forbidden region and thus is closer to the no-boundary proposal. Whether or not there are physically significant differences to the no-boundary wave function also at large volume remains to be studied.

This conclusion may seem bad if the no-boundary proposal does in fact not provide enough inflation. However, when we consider the evolution of a scalar or other matter field coupled to the discrete quantum geometry, there exist small-volume quantum modifications to the kinetic term of the matter Hamiltonian which automatically lead to inflation.Inflation  The modifications arise from the quantization of the inverse volume which removes the classical divergence In the situation of chaotic inflation where one couples an inflaton $`\varphi `$, the modified inflaton dynamics implies that $`\varphi `$ is driven up its potential at early stages (see Fig. 4). Before scales much larger than the Planck length are reached, the inflaton turns around and enters a slow-roll phase of chaotic inflation. Thus, the traditional considerations of probability distributions for the inflaton are not even needed here since there is a dynamical mechanism which results in large inflaton values. This occurs for a wide range of initial conditions with the initial value of the scalar field and its conjugate momentum being constrained by the uncertainty principle.InflationWMAP  The particular behavior may even have observable consequences since the slow-roll condition is violated around the turning point, i.e. at early stages of the slow-roll phase translating into large current scales.InflationWMAP

## IV Conclusions

The overall picture we have seen from loop quantum cosmology is very different from the traditional Wheeler–DeWitt quantization, and its structure is much richer. There are many phenomenological effects: A cut-off of classically diverging densities,InvScale ; Ambig  which is well-studied in flat models, implies inflationInflation  with possibly observable effectsInflationWMAP  when coupled to chaotic inflation. Also without coupling an inflaton, the universe would still expand with positive acceleration. Whether or not this is physically viable is undecided at this time since the usual perturbation theory becomes inapplicable.InflationWMAP

In the closed model with its classically forbidden region there are even more diverse phenomenological possibilities which have not yet been studied in detail. In this case, not just densities, or extrinsic curvature, experience a cut-off by quantum geometry effects, but also the intrinsic curvature can become suppressed at small volumes.Closed  There the closed model behaves similarly to the flat model with intrinsic curvature close to zero and a classically allowed region can be created near the singularity. It remains to be seen what physical consequences these effects have as the classically allowed regions lies deep within the quantum geometry regime at small volume.

Other open issues are concerned with details of the evolution through the classical singularity. Since this involves the dynamics deep inside the Planck regime, results are very sensitive to quantization choices. This may be used to find restrictions among possible choices, but there are also many different cases which have to be studied. The equation (III.1) used here implies a symmetric solution which leads to equal behavior at positive or negative $`\mu `$. But there is a more complicated versionIsoCosmo  which is closer to the constraint in the full theory.QSDI  The properties discussed here are insensitive to which one of those two forms for the constraint is used. But if we compare the wave function at positive and negative $`\mu `$ obtained from the more complicated constraint, then there will usually be significant small-scale oscillations at one side if they are constrained to be small at the other (see Fig. 5).DynIn  These oscillations arise because the more complicated constraint is of higher order in internal time than the simpler one (also $`\psi _{\mu \pm 8}`$ are involved). Being discrete, the evolution equation may contain coefficients which do not behave smoothly especially near the classical singularity. Thus passing through the singularity a smooth wave function can have the Planck scale frequency modes excited whose amplitudes can be very large in some models.ScalarLorentz  Particular quantizations choices can be favored if they provide an evolution equation with smooth coefficients, which may be exploited to reduce the freedom in the full quantization.

From the continuum picture point of view, the freedom for solutions to the discrete equation is reduced compared to the Wheeler–DeWitt equation. Since one has to require the existence of a continuum limit (or at least a good continuum approximationFundamentalDisc ) for the comparison, equation (III.1) is reduced to second order (even for the more complicated version mentioned in the preceding paragraph). The dynamical initial condition then leads to a reduction in the number of free parameters as seen before. From the discrete point of view, however, there are many more solutions (infinitely many since $`\mu `$ is continuous) with Planck scale oscillations. Their role and physical interpretation are open issues, and it is expected that an understanding requires the so far unknown physical inner product and the consideration of quantum observables. Tentative ideas in this direction can be found in \[cosmoIV, ; Golam, \], and further investigations are in progress.

There are promising possibilities for future developments. It is possible to extend the results to less symmetric models, which has already been completed for anisotropic but still homogeneous onesHomCosmo ; Spin  resulting in new physical effects.NonChaos  For cosmological considerations it is also important to add inhomogeneities at least as perturbations or in midisuperspace models. This will be done in the future in order to complete the picture provided by loop quantum cosmology, but even before these technical developments much remains to be studied, as indicated above, even in isotropic models.

M. B. is grateful to C. Kiefer for an invitation to a talk at the Xth Marcel Grossmann meeting, July 20–26, 2003, Rio de Janeiro, on which this paper is based. This work was supported in part by NSF grant PHY00-90091 and the Eberly research funds of Penn State.


text:
# Dynamical mean-field theory of Hubbard-Holstein model at half-filling: Zero temperature metal-insulator and insulator-insulator transitions

## Abstract

We study the Hubbard-Holstein model, which includes both the electron-electron and electron-phonon interactions characterized by $`U`$ and $`g`$, respectively, employing the dynamical mean-field theory combined with Wilson’s numerical renormalization group technique. A zero temperature phase diagram of metal-insulator and insulator-insulator transitions at half-filling is mapped out which exhibits the interplay between $`U`$ and $`g`$. As $`U`$ ($`g`$) is increased, a metal to Mott-Hubbard insulator (bipolaron insulator) transition occurs, and the two insulating states are distinct and can not be adiabatically connected. The nature of and transitions between the three states are discussed.

The interaction-induced phenomena are a fundamental problem because of the new physics that emerges. For instance, the interaction-driven metal-insulator transition (MIT) has been studied to a great extent for systems which involve the electron-electron or electron-phonon interactions imada . In real condensed-matter systems, both interactions exist and it will be important to understand the interplay between them. The prototype model for describing such interplay is the Hubbard-Holstein (HH) model, which includes the onsite electron-electron and electron-phonon interactions characterized by $`U`$ and $`g`$, respectively, as given by Eq. (1) below. The Hubbard model ($`g=0`$) and Holstein model ($`U=0`$), which are special limits of the more general HH model, have been extensively studied using the dynamical mean-field theory (DMFT) in infinite dimensions metzner ; georges  in the context of the interaction-driven MIT. Both models exhibit MIT as $`U`$ or $`g`$ is increased above a critical value. Now, consider the MIT in the $`Ug`$ plane. Natural questions then arise as to whether and how the two insulating states, metallic states, and metal-insulator transitions of the Hubbard and Holstein models are different, and how they are affected when both $`U`$ and $`g`$ are present and compete with each other. These questions will be addressed here by studying the HH model within the DMFT combined with Wilson’s numerical renormalization group (NRG) bulla2 . The NRG technique krishna1 ; wilson1  is particularly powerful in that it is non-perturbative in nature so that it can cover the whole parameter space and that it can probe the extremely small energy scales like the narrow coherence peak and soft phonon mode which emerge as the metal-insulator transitions are approached.

We will present a zero temperature phase diagram of unbroken symmetry ground states with a focus on the nature of and transitions between the ground states as $`U`$ and $`g`$ are varied. The main results of the present work are: (i) The ground state is a metallic (M) state when both $`U`$ and $`g`$ are small, but is a bipolaron insulating (BPI) state when $`g`$ is large, and is a Mott-Hubbard insulating (MHI) state when $`U`$ is large (see Fig. 1). In the shaded region in Fig. 1, the insulating and metallic solutions coexist, and the ground state is the one with a lower energy. (ii) The phase transition between M and BPI occurs along the dashed line that cuts through the coexistence region connecting E and T of Fig. 1. The order parameter ($`<n1>`$, see below.) changes discontinuously at the transition and the M-BPI transition is first order (Fig. 2). (iii) The transition between M and MHI occurs along the outer boundary of the coexistence region denoted by the solid line connecting T and U2. It is second order as in the Hubbard model (Fig. 3). (iv) The two kinds of insulators (BPI and MHI) are distinct and can $`not`$ be adiabatically connected. The transition between them is first order, and occurs along $`U=2g^2/\omega _0`$ (Fig. 4. See below). These points will be discussed in more detail below.

The HH model is defined by

$``$ $`=`$ $`{\displaystyle \frac{t}{\sqrt{q}}}{\displaystyle \underset{<i,j>\sigma }{}}c_{i\sigma }^{}c_{j\sigma }+U{\displaystyle \underset{i}{}}n_in_i`$ (1)
$`+`$ $`\omega _0{\displaystyle \underset{i}{}}a_i^{}a_i+g{\displaystyle \underset{i}{}}\left(a_i^{}+a_i\right)\left(n_i1\right),`$
where $`q`$ is the nearest-neighbor coordination number, $`n_{i\sigma }=c_{i\sigma }^{}c_{i\sigma }`$ the electron density operator for spin $`\sigma `$ at site $`i`$, and $`n_i=_\sigma n_{i\sigma }`$. The electrons have the onsite Coulomb repulsion $`U`$, and are linearly coupled with the Einstein phonon of frequency $`\omega _0`$ with the onsite coupling constant $`g`$. The DMFT was employed to solve the HH model of Eq. (1). The DMFT maps a lattice model onto an effective single site impurity model imbedded in a bath which is determined self-consistently via iterations metzner ; georges . The effective impurity problem was solved by the NRG technique in the present study.

The NRG+DMFT has been successfully applied to the half-filled Hubbard model to study the MIT at both $`T=0`$ bulla2  and finite temperatures bulla3 , and to the half-filled Holstein model at $`T=0`$ meyer . For the Hubbard model, the ground state is insulating for $`U>U_{c2}`$ and metallic for $`U<U_{c1}`$. Between $`U_{c1}<U<U_{c2}`$, the metallic and insulating solutions coexist, and the metallic state has the lower energy and is the ground state. The MIT occurs at $`U_c=U_{c2}`$, and is 2nd order. For the Holstein model, the ground state is insulating for $`g>g_{c2}`$ and metallic for $`g<g_{c1}`$. Meyer et al. reported that $`g_{c2}g_{c1}`$ is reduced as the phonon frequency $`\omega _0`$ is decreased and $`g_{c1}=g_{c2}`$ for $`\omega _0=0.05W`$, where $`W=4t`$ is the bandwidth meyer . We generalized these lines of research and study the HH model at half-filling with the NRG+DMFT incorporating the improved method for calculating electron bulla-jpc98  and phonon spectral function jeon . We adopt the semi-circular density of states which is realized in the Bethe lattice in infinite dimensions.

Our central result is the Fig. 1, where we present a zero temperature phase diagram in the $`UE_{bp}`$ plane of unbroken symmetry ground states of the HH model at half-filling obtained from the NRG+DMFT calculations. $`E_{bp}=2g^2/\omega _0`$ is the bipolaron binding energy which is twice the polaron energy, $`E_{bp}=2E_p`$. We take the bandwidth $`W`$ as the unit of energy ($`W=1`$). We take $`\omega _0=0.05`$ as with Meyer et al. Three distinct states are found in the $`UE_{bp}`$ plane: metal, bipolaron insulator, and Mott-Hubbard insulator states. We will look at them in more detail below. Specifically, we will address (i) the nature of M, BPI, and MHI states, (ii) transition between M and BPI, (iii) transition between M and MHI, and (iv) transition between the BPI and MHI. The results obtained from the NRG+DMFT regarding these points were summarized in the introduction above.

(i) Nature of M, BPI, and MHI states. – The ground state is a metallic state when both $`U`$ and $`g`$ are small. The M state becomes more correlated as $`g`$ or $`U`$ is increased as reflected in the decreasing quasiparticle weight $`z=\left[1\mathrm{\Sigma }_1(\omega )/\omega \right]_{\omega 0}^1,`$ where the subscript 1 (2) refers to the real (imaginary) part. $`z`$ vanishes at the outer boundary of the coexistence region. The BPI state is insulating because the attractive interaction between the electrons mediated by the phonons binds two electrons into a bipolaron as Capone and Ciuchi reported capone1 . The MHI state, on the other hand, is insulating because of the repulsion $`U`$ between two electrons at the same site.

An important point for understanding the overall feature of the phase diagram is that the effective interaction between two electrons in the HH model is given, after integrating the phonons out, by

$`V_{eff}(\omega )=U+{\displaystyle \frac{2g^2\omega _0}{\omega ^2\omega _0^2}}.`$ (2)
The interactions between electrons due to $`g`$ and $`U`$ compete each other for $`\omega \omega _0`$, which yields $`V_{eff}(\omega =0)=0`$ along $`U=E_{bp}`$. The $`U=E_{bp}`$ line plays a special role in the phase diagram in that all three states meet at the triple point T lying on this line and the BPI-MHI transition occurs along it, as we will discuss below. Also, the critical $`g_c`$, which separates M and BPI along E and T, should increase as $`U`$ is increased because the attractive interaction due to $`g`$ has to overcome the repulsive interaction of $`U`$ to bind two electrons into a bipolaron. This expectation is indeed borne out as shown in Fig. 1. On the other hand, the $`U_c`$, which is the solid line connecting T and U2, depends weakly on $`g`$.

For the coexistence region with $`U>E_{bp}`$, the arguments which led to the conclusion of the MIT of the Hubbard model being a 2nd order at $`U=U_{c2}`$ moeller  still apply and the M is the ground state. On the other hand, in the coexistence region with $`E_{bp}>U`$, the $`<n_{}n_{}>_{BPI}`$ and $`<n_{}n_{}>_M`$, where $`<𝒪>_{M(I)}`$ stands for an expectation value of an operator $`𝒪`$ in the metallic (insulating) state, do not meet as shown in Fig. 2(a). This implies that the M-BPI transition is 1st order. The energy calculations in Fig. 2(a) show that the M and BPI solutions cross each other inside the coexistence region.

(ii) Transition between M and BPI. – The MIT of half-filled Holstein model in infinite dimensions was discussed previously by Meyer et al. meyer  and by Benedetti and Zeyher benedetti . The BPI state is insulating because electrons are bound in bipolarons. Because we set the chemical potential $`\mu =U/2`$, which prefers the electron density $`<n>=1`$, the bipolaron formation may be accommodated by reconstructuring the system into a phase separated state or a charge ordered state where the doubly occupied and empty sites alternate in the real space. This possibility is not allowed in the present DMFT work because the unit cell consists of a single site. Instead, the bipolaron instability shows up as degenerate ground states: In the NRG+DMFT calculations, the ground state is 8 fold degenerate in the BPI regime, with a one-to-one correspondence between the two sets of $`<n>11`$ and $`<n>11`$, which on averaging satisfies $`<n>=1`$. The order parameters for the M-BPI may be taken as

$`C=<n1>.`$ (3)
Now, consider the cut 2b in Fig. 1 along which we increase $`g`$ with $`U=1.0`$. Let $`g_{c1}`$ and $`g_{c2}`$ be, respectively, the lower and upper boundaries of the shaded coexistence region. As $`gg_{c2}`$, the quasiparticle weight $`z`$ and the renormalized phonon frequency $`\mathrm{\Omega }`$ approach 0 continuously, but, the order parameter $`C`$ increases from 0 to a finite value discontinuously, as shown in Fig. 2(b) and (d). The softening phonon mode is a manifestation of a lattice instability just like structural phase transitions. A stability is restored by a condensation of the unstable mode. It results in a non-zero expectation value of the phonon operator ($`<a>0`$), which may be taken as an order parameter for the M-BPI transition. The phonon hardens back to the bare mode as $`g`$ is increased above $`g_{c2}`$ because the screening is not effective in an insulating state. The $`a`$ and $`(n1)`$ are linearly coupled in the Hamiltonian of Eq. (1), and a zero/nonzero $`<a>`$ is expected to lead to a zero/nonzero $`C`$. Therefore, an equally good order parameter is the $`C`$ of Eq. (3).

In Fig. 2(b), we show $`z`$, $`C`$, and $`\frac{\omega _0}{g}<a>`$ with increasing $`g`$. Note that $`<n1>\frac{\omega _0}{g}<a>`$ as expected from a simple mean-field theory. In Fig. 2(c), we plot the renormalized electron spectral function, $`A(\omega )=\frac{1}{\pi }Imc,c^{}_\omega `$, to show the MIT as $`g`$ is increased along the cut 2b. In Fig. 2(d), we plot the renormalized phonon spectral function, $`B(\omega )=\frac{1}{\pi }Ima+a^{},a+a^{}_\omega `$, calculated as $`g`$ is increased using the method in jeon . Note that at $`g=0.1664`$ which is close to $`g_{c2}=0.1665`$, the phonon mode splits into two components as reported by Jeon $`etal.`$ for the Anderson-Holstein model: one component develops into the soft mode and the other hardens back to the bare mode jeon . But, in the coexistence region where $`g_{c1}`$ and $`g_{c2}`$ are not close, the 1st order MIT from M to BPI preempts the emergence of the soft mode as $`g`$ is increased because the MIT occurs at $`g_c`$ ($`g_{c1}<g_c<g_{c2}`$).

The M-BPI transition may be understood from the Ginzburg-Landau theory. The free energy may be written in terms of the the local phonon coordinate $`\xi =\sqrt{\frac{\mathrm{}}{2M\omega _0}}(a^{}+a)`$ as

$`F={\displaystyle \frac{1}{2}}\alpha \xi ^2+{\displaystyle \frac{1}{4}}\beta \xi ^4+{\displaystyle \frac{1}{6}}\gamma \xi ^6,`$ (4)
with $`\gamma >0`$. The phase transition can be 1st order for $`\beta <0`$. In the context of the M-BPI transition, the $`g_{c1}<g<g_{c2}`$ regime corresponds to $`\beta ^2/4\gamma >\alpha >0`$, for which three local minima exit as a function of $`\xi `$. The $`g_{c1}`$, $`g_c`$, and $`g_{c2}`$ correspond, respectively, to $`\alpha =\beta ^2/4\gamma `$, $`3\beta ^2/16\gamma `$, and 0. Although the $`\alpha `$ and $`\beta `$ are complicated functions of $`U`$ and $`g`$, the DMFT phase diagram of Fig. 1 indicates that $`\alpha `$ and $`\beta `$ are mainly determined by, respectively, $`g`$ and $`U`$. This observation suggests an interesting possibility that the M-BPI transition may become 2nd order for small values of $`U`$ where $`\beta 0`$, and a quantum tricritical point exists between E and T. This possibility will be explored further in a subsequent study.

(iii) Transition between M and MHI. – As the MIT from M to BPI is driven by growing charge susceptibility, the one from M to MHI is driven by growing spin susceptibility. In Fig. 3(a), we plot $`z`$ as $`U`$ is increased, and $`<n_{}n_{}>`$ as $`U`$ is increased and decreased for a fixed $`g=0.1`$. $`z`$ goes to 0 continuously, and $`<n_{}n_{}>_M=<n_{}n_{}>_{MHI}`$ as $`UU_{c2}`$. The M-MHI transition at $`U_c=U_{c2}`$ is, therefore, 2nd order. The phonons are hardly renormalized (not shown) because the charge degree of freedom is not soft unlike the MIT between M and BPI. In Fig. 3(b), $`A(\omega )`$ is plotted to show the MIT as $`U`$ is increased.

The $`U_c`$ does depend on $`g`$, although weakly, as can be seen from the Fig. 1, which implies that the M-MHI transition is coupled with the phonons. This raises an interesting possibility that some discrepancies between the experimental observations on V2O3 systems and DMFT results of the Hubbard model may be resolved in terms of the electron-phonon coupling limelette .

(iv) Transition between BPI and MHI. – Now, we turn to an interesting question about the two kinds of insulators: Are BPI and MHI different, or, put differently, can the two insulating states be adiabatically connected? Along the line $`U=E_{bp}`$, the effective interaction in the static limit vanishes. One might then naively expect that a metallic state is the ground state along this line all the way to $`U=E_{bp}\mathrm{}`$. One can show, however, that for $`U=E_{bp}\mathrm{}`$ the ground state is an insulating state using the exact atomic limit Green’s function jeon ; hewson , as found from the detailed NGR+DMFT calculations.

In Fig. 4, we show $`C`$ as $`g`$ is varied for $`U=2.5`$. The $`C`$ changes discontinuously exactly at $`E_{bp}=U`$ and does not exhibit any other non-analytic behavior as $`g`$ is varied. We checked another point on the line $`E_{bp}=U`$ by varying $`U`$ with a fixed $`g=0.3`$. The transition occurs exactly at $`E_{bp}=U`$ (not shown). This confirms that the two insulating states are distinct, and the phase transition between them occurs along the line of $`U=E_{bp}`$ and is 1st order. The point T in Fig. 1 is a triple point where two 1st order and one 2nd order phase transition lines meet.

We would like to thank Tae-Suk Kim, Yunkyu Bang, Ralf Bulla, Alex Hewson, and Dietrich Meyer for helpful comments and discussions. This work was supported by the Korea Science & Engineering Foundation (KOSEF) through grant No. R01-1999-000-00031-0 and Center for Nanotubes and Nanoscale Composites (CNNC), and by the Ministry of Education through BK21 SNU-SKKU Physics program. JHH and HCL acknowledge the support from Center for Strongly Correlated Materials Research (CSCMR).


text:
# Eccentricity fluctuations and its possible effect on elliptic flow measurements

## Abstract

The elliptic flow measured at RHIC has been interpreted as a signature for strong partonic interactions early in the collision and as an indication of a well developed quark-gluon plasma phase. The measured values of elliptic flow, using methods based on multi-particle correlations, are affected by fluctuations in the magnitude of the elliptic flow. In this Letter, using a Monte Carlo Glauber calculation, we estimate what the possible effect of spatial eccentricity fluctuations is on the determination of elliptic flow.

In non-central heavy-ion collisions, the initial spatial anisotropy due to the geometry of the overlap region and the pressure developed early in the collision generate an observable azimuthal momentum-space anisotropy. The particle yields produced in heavy-ion collisions can be characterized by VoloshinFourier :

$$\frac{\mathrm{d}^3N}{\mathrm{d}p_t^2\mathrm{d}\varphi \mathrm{d}y}=\frac{\mathrm{d}^2N}{2\pi \mathrm{d}p_t^2\mathrm{d}y}[1+2\underset{n}{}v_n\mathrm{cos}(n(\varphi \mathrm{\Psi }_R))],$$
(1)
where $`p_t`$ is the transverse momentum of the particle, $`\varphi `$ is its azimuthal angle, $`y`$ is the rapidity and $`\mathrm{\Psi }_R`$ the reaction plane angle, see fig 1. The second coefficient, $`v_2`$, of this Fourier series is called elliptic flow.

Elliptic flow as a signature of hydrodynamic behavior of nuclear matter produced in high energy nuclear collisions has been proposed by Ollitrault in 1992 Ollitrault92 . After that it has attracted increasing attention from both experimentalists and theorists OverviewFlowTheory  and has been measured at AGS E877 ; E895 , SPS NA49 ; CERES  and RHIC STAR ; PHOBOS ; PHENIX  energies. It is thought that elliptic flow reflects the amount of interactions between the constituents at an early time in the evolution of the produced system SorgeEarly . Therefore it is sensitive to the equation of state of the produced system when this system might be in the quark-gluon plasma phase.

Since the reaction plane is not known experimentally, the elliptic flow is calculated using azimuthal angular correlations between the observed particles PoskanzerVoloshinMethods . In the case of two particle correlations the measurement is proportional to $`v_2^2`$. The reported elliptic flow values are therefore obtained as $`\sqrt{v_2^2}`$ after averaging over events.

Because elliptic flow is a collective effect, it is a correlation of all the particles with the reaction plane. This can be exploited experimentally by using multiple particle correlations to calculate $`v_2`$. To calculate these correlations a convenient mathematical approach is to use cumulants. This method, proposed in OllitraultCumulants , has the additional advantage that it allows to subtract the so called non-flow effects from $`v_2`$. Non-flow effects are correlations between the particles not related to the reaction plane. Such effects include, but are not limited to, resonance decays, (mini)jet fragmentation and Bose-Einstein correlations. The cumulant method uses multi-particle correlations which introduce higher powers of $`v_2`$. The corresponding equations for calculating $`v_2`$ in the cumulant method for two, four and six particle azimuthal correlations are given by:

$`(v_2\{2\})^2`$ $`=`$ $`v_2^2`$
$`(v_2\{4\})^4`$ $`=`$ $`2v_2^2^2v_2^4`$
$`(v_2\{6\})^6`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left(v_2^69v_2^4v_2^2+12v_2^2^3\right),`$ (2)
where $`v_2\{2\}`$ is calculated using two-particle azimuthal correlations, $`v_2\{4\}`$ using a mix of two and four-particle azimuthal correlations and $`v_2\{6\}`$ uses two, four and six-particle azimuthal correlations.

However, due to event-by-event fluctuations in the elliptic flow for instance, the event averaged $`v_2^nv_2^n`$ for $`n2`$ STAR\_flowPRC . Therefore, comparing the experimental values of $`v_2`$ with model calculations which do not correctly include these fluctuations might not be a priory justified. Furthermore, there are effects due to the finite width of the centrality bins where also $`v_2^nv_2^n`$. The STAR collaboration STAR\_flowPRC  found a negligible bias on the extracted $`v_2`$ due to these binning effects, except in the most central bin.

To estimate the possible effect of the fluctuations on the measured elliptic flow values we calculate the fluctuations in the initial spatial anisotropy of the created system using a Monte Carlo Glauber model (MCG). This anisotropy, which generates the elliptic flow, is given by SorgeEpsilon ; Heiselberg :

$$ϵ\frac{y_i^2x_i^2}{y_i^2+x_i^2},$$
(3)
where $`x_i^{}`$ and $`y_i^{}`$ are the coordinates of the constituents in the plane perpendicular to the beam and $`x^{}`$ is in the reaction plane (see Fig. 1). Due to the relation between the initial spatial anisotropy and the elliptic flow, $`v_2`$ $`ϵ`$ Heiselberg ; SorgeEpsilon ; PoskanzerVoloshinEpsilon , fluctuations in $`ϵ`$ will lead to fluctuations in $`v_2`$. It should be noted that this is only one specific example of fluctuations which could contribute to the measured value of $`v_2`$.

The MCG approach allows for an event-by-event calculation of $`ϵ`$ and therefore the determination of the corresponding higher order moments of the event averaged distribution of $`ϵ`$. For details on the MCG model and the parameters used see Ref. STAR\_pionGlauber .

In Fig. 2a, the eccentricity calculated using the MCG is plotted versus the impact parameter ($`b`$) of an Au+Au collision. This figure shows that the fluctuations in $`ϵ`$ for the most peripheral collisions are large. For the most central collisions $`ϵ`$ can be both positive and negative. This leads to an obvious bias when calculating $`ϵ`$ from $`ϵ^2^{1/2}`$. Figure 2b shows the calculated $`ϵ`$ and the corresponding $`ϵ^n^{1/n}`$. It is clear that there is a bias over the whole centrality range, however this bias is the largest for the most central and most peripheral collisions.

Using $`v_2ϵ`$, we replace in Eq. 2 $`v_2`$ by $`ϵ`$, which allows us to calculate $`ϵ\{2\}`$, $`ϵ\{4\}`$ and $`ϵ\{6\}`$. If the right-hand side of Eq. 2 becomes negative, we take the $`n^{\mathrm{th}}`$ root of the absolute value and multiply this by $`1`$. Figure 3 shows the ratio of $`ϵ\{m\}`$ and $`ϵ`$, where $`m=2,4,6`$, versus the collision centrality in terms of cross-section. In the standard MCG approach for nucleons with a cross-section of 42 mb the event by event fluctuations lead to a ratio $`ϵ\{2\}/ϵ`$ which is always larger than 1. This shows that the experimental determination of $`v_2`$, using two-particle azimuthal correlations in the case of event-by-event fluctuations in $`ϵ`$, leads to an overestimation of the true elliptic flow. This is particularly true for the most central and most peripheral events. The inset in Fig. 3 shows the same calculated values on a expanded scale. Even for the centrality region of 20-60% the two-particle correlation method overestimates the true elliptic flow by about 10%. In the same figure the results for the higher order cumulant ratios, $`ϵ\{4\}/ϵ`$ and $`ϵ\{6\}/ϵ`$, are shown. The higher order cumulants in the case of event by event fluctuations underestimate the true elliptic flow in the centrality range of 0-80%. Above 80% of the cross section also the higher order cumulants overestimate the elliptic flow. In the mid-central region (20-60%) the higher order cumulants underestimate the true elliptic flow by about 10%. If the fluctuations are small it can be shown analitically that the true value of $`v_2`$ is given by $`v_2\{2\}+v_2\{4\}/2`$ Ollipriv .

The most central 10% of the cross section show an interesting behavior for the different cumulants. The two particle correlations always lead to a real value of $`ϵ\{2\}`$, however the higher order cumulants can become complex because of the combination of two, four and higher multi-particle contributions. This is the case for $`ϵ\{4\}`$ and $`ϵ\{6\}`$ in the centrality bin 5-10%. However, in the most central bin, 0-5%, $`ϵ\{6\}`$ becomes real again. Complex values of $`v_2`$ are usually not reported by experiments.

To investigate the sensitivity to the magnitude of the fluctuations in the eccentricity we followed the work of Eremin and Voloshin VoloshinQuarkPart  and calculate $`ϵ\{m\}`$ using valence quarks as constituents in the MCG. The cross section used for the quarks was set to 6 mb and the number of constituent quarks in a gold nucleus to 591. All the other parameters were kept the same as in the nucleon MCG. Due to the larger number of interacting constituents the fluctuations in central collisions are reduced. Indeed, the $`ϵ`$ fluctuations for the 0-80% of the cross-section are reduced as shown by the lines in Fig. 3.

The ratios shown in Fig. 3 cannot be directly compared to the elliptic flow measurements since $`v_2`$ is not an experimental observable. Instead, we calculate the ratio $`ϵ\{m\}/ϵ\{n\}`$. Because of the assumed proportionality between $`ϵ`$ and $`v_2`$, which in principle could depend on centrality STAR\_flowPRC , this ratio can be compared with the measured $`v_2\{m\}/v_2\{n\}`$.

Figure 4a shows the ratio $`ϵ\{4\}/ϵ\{2\}`$ for the quark and nucleon MCG calculations. In both cases this ratio is smaller than unity over the whole centrality range, with the largest suppression in the case of the nucleon MCG. The elliptic flow cumulant measurements from the STAR collaboration STAR\_flowPRC  are also shown in Fig. 4. The experimental results are in between the calculated values, and are closer to the nucleon (quark) MCG results for peripheral (central) collisions.

If the observed difference between $`v_2\{2\}`$ and $`v_2\{4\}`$ would be due to non-flow effects, it is proposed in Ref. OllitraultCumulants  that the squared difference is given by

$$v_2^2\{2\}v_2^2\{4\}=\frac{g_2}{N}$$
(4)
where $`N`$ is the multiplicity and $`g_2/N`$ is the non-flow contribution. NA49 NA49flow  has indeed observed such a constant behavior.

This squared difference, as measured by the STAR collaboration STAR\_flowPRC , is shown in Fig. 4b for the cumulant measurements as a function of centrality. Instead of scaling with the multiplicity as done by NA49 we scaled the STAR results using the corresponding number of wounded nucleons ($`N_W`$), which is proportional to the multiplicity. These measurements are compared with the values obtained from the quark and nucleon MCG (histograms in Fig. 4b). However, to make this comparison the proportionality between $`v_2`$ and $`ϵ`$ was assumed to be 0.16, independent of centrality. It is seen from the figure that the centrality dependence of the data is better described by the nucleon MCG than by a constant value as expected for non-flow contributions.

In summary, experimental measurements of elliptic flow ($`v_2`$) might be affected by fluctuations. In this study we investigate the effect on the determination of $`v_2`$ of fluctuations in the spatial anisotropy of the overlap region ($`ϵ`$) in heavy-ion collisions using a Monte Carlo Glauber simulation. For this study we assume a proportional relationship between $`ϵ`$ and $`v_2`$.

When elliptic flow is measured from two-particle correlations ($`v_2\{2\}`$) we find that fluctuations in $`ϵ`$ lead to an overestimation of its value. On the other hand a cumulant analysis of order four ($`v_2\{4\}`$) and six ($`v_2\{6\}`$) leads to an underestimation of the true value.

Experimental results on the ratios $`v_2\{m\}/v_2\{n\}`$ can be compared directly to the model and are shown to be in approximate agreement. Traditionally, a ratio different from unity is interpreted as being due to non-flow effects. However, this study shows that this difference can also be explained by event-by-event fluctuations in the elliptic flow. In reality it is likely that both non-flow and fluctuations affect the measured elliptic flow. It remains to be determined which effect dominates.

Experimentally, this can be addressed by also reporting the ratios $`v_2\{m\}/v_2\{n\}`$ obtained with the cumulant analysis as a function of centrality when they become complex because this is a unique feature of fluctuations. Event-by-event fluctuations in elliptic flow have been suggested as an interesting physics observable otherFluct , however to address this these more “trivial” causes of fluctuations have to be understood first.
Acknowledgments:

We are grateful to Sergei Voloshin, Aihong Tang, Art Poskanzer, Jean-Yves Ollitrault, Peter Kolb, Pasi Huovinen and Michiel Botje for useful discussions and suggestions.