text:
# Small Fermi surface in the one-dimensional Kondo lattice model

(March 18, 2024)

## Abstract

We study the one-dimensional Kondo lattice model through the density matrix renormalization group (DMRG). Our results for the spin correlation function indicate the presence of a small Fermi surface in large portions of the phase diagram, in contrast to some previous studies that used the same technique. We argue that the discrepancy is due to the open boundary conditions, which introduce strong charge perturbations that strongly affect the spin Friedel oscillations.

Several uncertainties still exist in our understanding of the physics of heavy fermion materials.A. C. Hewson (1993) The importance of solving these uncertainties has become even more pressing as we attempt to understand the anomalous behavior observed in the vicinity of the clean antiferromagnetic quantum critical point.P. Coleman et al. (2001) The two major scenarios take radically different points of view. In the first one, conduction electrons are assumed to acquire their peculiar dynamics through an essentially perturbative coupling to the slow critical modes of the antiferromagnetic background.J. A. Hertz (1976); A. J. Millis (1993) Alternatively, the starting point is taken to be the strong coupling of the conduction electrons and the localized spins to form singlets, as in the single impurity Kondo problem. The nature of the quantum critical point is then linked to a non-trivial competition between the local dynamics and the long wavelength antiferromagnetic fluctuations.Q. Si et al. (2001)

One of the key assumptions of the second approach is the presence of a large Fermi surface (FS), namely one whose volume is given by including the localized spins in the count

$$\frac{2V_{FS}^{}}{\left(2\pi \right)^D}=n+1.$$

Behind this assumption lies a deep connectionR. M. Martin (1982) between the Friedel sum rule of the single impurity Kondo problemJ. Friedel (1952) and Luttinger’s theorem for the FS volume of a system of interacting fermions.J. M. Luttinger (1960) The inclusion of the localized electron in the count is natural within an Anderson lattice model description at weak coupling but becomes doubtful at strong coupling where the Kondo lattice Hamiltonian is the effective low-energy theory. However, topological arguments have been used, in oneM. Yamanaka et al. (1997) as well as in higherM. Oshikawa (2000) dimensions, to show that indeed neutral gapless excitations must exist at a $`𝐤`$-vector spanning a large FS.111We adhere to the somewhat inappropriate yet common usage of the term Fermi “surface” in one dimension.  Furthermore, numerical studies of the one-dimensional Kondo lattice model have also pointed to the presence of a large FS.N. Shibata et al. (1996, 1997) In this paper, we show that a more careful analysis of the numerical evidence casts serious doubts on these conclusions and leaves open the question of the size of the Fermi surface in heavy fermion systems.

We considered the one-dimensional $`S=\frac{1}{2}`$ Kondo lattice Hamiltonian with $`L`$ sites

$$H=t\underset{j=1,\sigma }{\overset{L1}{}}c_{j,\sigma }^{}c_{j+1,\sigma }^{}+J\underset{j=1}{\overset{L}{}}𝐒_j𝐬_j$$

where $`c_{j\sigma }`$ annihilates a conduction electron in site $`j`$ with spin projection $`\sigma `$, $`𝐒_j`$ is a localized spin $`\frac{1}{2}`$ operator and $`𝐬_j=\frac{1}{2}_{\alpha \beta }c_{j,\alpha }^{}𝝈_{\alpha \beta }c_{j,\beta }^{}`$ is the conduction electron spin density operator. $`J>0`$ is the Kondo coupling constant between the conduction electrons and the local moments and the hopping amplitude $`t`$ is set to unity to fix the energy scale. We treated the model with the density matrix renormalization group (DMRG) techniqueS. R. White (1992, 1993) with open boundary conditions. We used the finite-size algorithm for sizes up to $`L=120`$ keeping up to $`m=400`$ states per block. The discarded weight was typically about $`10^510^8`$ in the final sweep.

There is compelling evidenceH. Tsunetsugu et al. (1997) that the one-dimensional Kondo lattice model away from half-filling is a Tomonaga-Luttinger (TL) liquid.J. Voit (1994) TL liquids with periodic boundary conditions, have charge and spin correlation functions given asymptotically by

$`\delta n\left(0\right)\delta n\left(x\right)`$ $`=`$ $`{\displaystyle \frac{K_\rho }{\left(\pi x\right)^2}}+A_1{\displaystyle \frac{\mathrm{cos}\left(2k_Fx\right)}{x^{K_\rho +1}}}`$ (1)
$`+`$ $`A_2{\displaystyle \frac{\mathrm{cos}\left(4k_Fx\right)}{x^{4K_\rho }}},`$
$`𝐒\left(0\right)𝐒\left(x\right)`$ $`=`$ $`{\displaystyle \frac{1}{\left(\pi x\right)^2}}+B_1{\displaystyle \frac{\mathrm{cos}\left(2k_Fx\right)}{x^{K_\rho +1}}},`$ (2)
where $`K_\rho `$ is a non-universal correlation exponent and $`k_F`$ is the Fermi momentum. Besides, local charge and spin perturbations, such as introduced by impurities or boundaries, lead to corresponding Friedel oscillations, which in the case of a TL liquid take the formN. Shibata et al. (1996, 1997); S. Eggert and I. Affleck (1995); R. Egger and H. Grabert (1995); J. Voit et al. (2000)

$`\delta n\left(x\right)`$ $`=`$ $`C_1{\displaystyle \frac{\mathrm{cos}\left(2k_Fx\right)}{x^{\left(K_\rho +1\right)/2}}}+C_2{\displaystyle \frac{\mathrm{cos}\left(4k_Fx\right)}{x^{2K_\rho }}},`$
$`\delta S_z\left(x\right)`$ $`=`$ $`D_1{\displaystyle \frac{\mathrm{cos}\left(2k_Fx\right)}{x^{K_\rho }}}.`$
We first show our results for the total spin correlation function $`c(j,k)=𝐒_j^T𝐒_k^T`$, where $`𝐒_j^T=𝐒_j+𝐬_j`$. Since we use open boundary conditions, translational invariance is lost and $`c(j,k)`$ depends on both $`j`$ and $`k`$. We present results obtained by two different methods. In the first one (dashed lines in Fig. 1), $`c\left(l=\left|jk\right|\right)`$ was obtained by taking $`j`$ and $`k`$ at the same distance (within a lattice spacing) from the center of the chain. We call this the central value of $`c\left(l\right)`$. In the second one (solid lines in Fig. 1), we averaged over all pairs of sites separated by the distance $`l=\left|jk\right|`$. We will call it the average value of $`c\left(l\right)`$. If the boundary perturbation has a negligible effect on the spin-spin correlation function, then the two methods should yield similar results and we can have confidence that we have obtained the bulk value of $`c\left(l\right)`$. This is indeed what is observed for the density $`n=\frac{2}{5}`$ with $`L=60`$, as seen in Fig. 1. The paramagnetic (PM) and ferromagnetic (FM) phasesH. Tsunetsugu et al. (1997) can be easily discerned from the long distance behavior of $`c\left(l\right)`$: for $`J<J_c`$, the envelope of $`c\left(l\right)`$ tends to zero (main plots in Fig. 1) and for $`J>J_c`$ it tends to the magnetization squared (inset of Fig. 1). The critical value is $`J_c1.2`$ for $`n=\frac{2}{5}`$. The ferromagnetism was also checked directly from the total spin of the ground state. The PM phase exhibits well developed Friedel oscillations due to the open boundaries.

In Fig. 2(a), we show the spin structure factor $`S\left(q\right)=\frac{1}{L}_{j,k}e^{iq\left(jk\right)}c(j,k)`$ corresponding to Fig. 1. Because of the weak influence of the boundaries, this is very close to the Fourier transform of $`c\left(l\right)`$. Whereas $`S\left(q\right)`$ is maximum at $`q=0`$ in the FM phase (with $`LS\left(0\right)=S^T\left(S^T+1\right)`$), in the PM phase $`S\left(q\right)`$ is peaked at $`q_s=\pi n`$ at $`n=\frac{2}{5}`$.222The residual value of $`S\left(q=0\right)`$ in the PM phase is due to the DMRG truncation at $`m=400`$. We have checked that it decays exponentially to zero as $`m`$ is increased.  This result does not seem to be due to finite size effects. Indeed, the structure factor peak gets sharper and more pronounced as one goes from $`L=60`$ to $`L=120.`$ This is shown in Fig. 2(b), for the two densities $`n=\frac{2}{5}`$ and $`n=\frac{4}{5}`$ at $`J=0.35.`$ This is *strong evidence for a small Fermi surface* with $`2k_F=q_s=\pi n`$, in sharp contrast to what was previously reported.N. Shibata et al. (1996, 1997) Note a very small feature at wave vector $`2\pi n`$, which, however, does not increase with the system size and is below the accuracy of the DMRG. We also calculated the spin structure factor at several other density values, as shown in Fig. 3(a). As finite size effects appear to be negligible, we have kept to smaller chain sizes ($`L=40`$ or $`30`$). In all cases, there was good agreement between central and average values of $`c\left(l\right)`$. All the results point to the presence of a small FS. In order to understand the origin of this discrepancy we re-examined the parameter ranges studied in Refs. N. Shibata et al., 1996 and N. Shibata et al., 1997. As we will see, their results occur at large values of $`J`$ ($`J1`$), where strong boundary charge perturbations mask the true bulk behavior of spin correlations. This is diagnosed by very different values of the average and the central $`c\left(l\right)`$. On the other hand, when $`J1`$ as in all cases seen above, the boundaries induce only a weak charge disturbance. As a result, central and average $`c\left(l\right)`$ are the same, and the spin correlation function oscillates at $`2k_F=q_s=\pi n`$.

In Fig. 3(b), we present the local density $`n_j`$ and the squared $`z`$-component of the total spin $`\left(S_{jz}^T\right)^2`$ versus distance for the density $`n=\frac{6}{7}`$ and $`J=2.5`$. The charge Friedel oscillations induced by the open boundaries are the same as found previouslyN. Shibata et al. (1996, 1997) and the squared $`z`$-component is anti-correlated with the charge. This strong charge disturbance is seen at other densities when $`J1`$ but is suppressed as $`J`$ is decreased. This can be seen in Fig. 4, where we show the magnitude of the Fourier transform $`N\left(q\right)=\frac{1}{L}_je^{iqj}\left(n_jn\right)`$ versus momentum for several densities and coupling constants. Furthermore, when $`J`$ is decreased, the charge Friedel oscillation peak moves from $`q_c=2\pi \left(1n\right)`$ to $`q_c=2k_F=\pi n`$. We note that a peak at $`q_c=2\pi \left(1n\right)`$ (mod $`2\pi `$) cannot distinguish between a small FS, with $`4k_F=2\pi n`$, and a large one, with $`4k_F^{}=2\pi \left(1+n\right)`$ (we denote a large Fermi vector by a star). For this reason, we cannot use the charge structure factor to determine the size of the Fermi surface. Even when the peaks are not sharp, the oscillations are still quite well defined, as shown in Fig. 5(a) for $`n=\frac{1}{2}`$ (compare the scales of Fig. 5(a) and Fig. 3(b)).

The presence of strong boundary charge disturbances leads to different behaviors of the central and the average spin correlation functions. This connection is made apparent in Figs. 5(b) and 6(a), which show average and central $`c\left(l\right)`$ for $`n=\frac{1}{2}`$, $`J=0.5`$ and $`n=\frac{4}{5}`$, $`J=1.5`$, respectively. At $`n=\frac{1}{2}`$, $`J=0.5`$, charge oscillations are small and the average and central $`c\left(l\right)`$ are almost the same. On the other hand, the large charge peak that occurs at $`n=\frac{4}{5}`$, $`J=1.5`$ (Fig. 4) leads to quite different values of the average and central $`c\left(l\right)`$. This is also reflected in the spin structure factor, which shows only broad ill-defined features, as plotted in Fig. 6(b). This is the key to understanding the discrepancy between our results and the ones of Shibata *et al.*N. Shibata et al. (1996, 1997)

In order to observe spin oscillations, Shibata et. al applied a small local field at the chain ends and measured $`S_{jz}^T`$. We obtained similar spin oscillations for $`\left(S_{jz}^T\right)^2`$, but with half the period, *without any field*. The origin of the latter structure is clear: local spin fluctuations are determined by the charge Friedel oscillations (see Fig. 3(b)). Furthermore, this also shows that the system is rendered more polarizable by the boundary perturbation at the peaks of $`\left(S_{jz}^T\right)^2`$, accounting for the oscillations of $`S_{jz}^T`$ when a magnetic field is applied at the endpoints. Thus, the response of the system under the application of boundary magnetic fields, as was done in Refs. N. Shibata et al., 1996 and N. Shibata et al., 1997, is strongly perturbed by the presence of open boundaries and cannot by itself be used to determine the size of the FS. In that case, the spin oscillations are inextricably linked to the charge ones by $`q_s=q_c/2`$. Additional confidence in this picture can be gained by the inspection of Figs. 2 and 4 of Ref. N. Shibata et al., 1996, where an increase of the charge peak is accompanied by an increase of the spin peak. Spin Friedel oscillations should ideally be studied by applying a small local magnetic perturbation *in the absence of any charge perturbation*. We conclude that, based on the available evidence, it is impossible to determine whether the system has a large or a small Fermi surface for $`J1`$.

We can envisage three alternatives to try to overcome this difficulty. The first one would be to use larger boundary fields so that the spin perturbation can surpass the charge one. However, it is not clear that this regime is attainable without drastically changing the nature of the ground state. A second way would be to work with periodic boundary conditions and a magnetic field applied at one site only, thus eliminating boundary charge perturbations while keeping spin ones. However, this is not particularly appealing because the DMRG loses accuracy with periodic boundary conditions. Finally, by going to larger systems one can have access to the true bulk behavior. Without these further studies, the size of the Fermi surface at $`J1`$ remains undetermined.

We have also calculated the spin correlations between localized spins only and conduction electrons only. In Fig. 7(a) and (b) we present these correlations together with the total spin-spin correlations for densities $`n=1`$ and $`n=\frac{1}{4}`$, respectively. In Fig. 7(b), the correlations between localized spins are not shown since they differ very little from the total spin ones. As expected all three functions have the same period, differing only in amplitude. It is also clear that the conduction electron contribution increases with the density. Note also that the spin correlations at half filling (Fig. 7(a)) decay much faster than at other fillings, due the presence of a spin gap.H. Tsunetsugu et al. (1992); C. C. Yu and S. R. White (1993)

Let us now discuss these results. The theorem of Ref. M. Yamanaka et al., 1997 guarantees that, provided there is either unbroken time reversal or parity symmetries in the ground state, the one-dimensional Kondo lattice has gapless excitations at $`2k_F^{}=\pi \left(1n\right).`$ This should be valid within the PM phase. The conventional Luttinger liquid phenomenology then tells us that these are collective spin and/or charge bosonic excitations with a linear dispersion with respect to deviations from this wave vector. They lead to the characteristic oscillatory behavior of Eqs. (1) and (2) and the corresponding peaks in the spin and charge structure factors. Of course, the theorem does not forbid the appearance of gapless excitations at *other* wave vectors such as $`2k_F=\pi n`$. Together with our results, this would seem to indicate that the spectral weight at $`2k_F^{}`$ is rather small in the range $`J1`$, most of it being concentrated at $`2k_F`$.

The conventional wisdom about heavy fermion systems is based on the competition between the local Kondo effect and the inter-site Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction. Although this dichotomy is controversial in one dimension, it is tempting to use it as a general guide. The size of the compensating Kondo cloud around a single localized moment has been argued to be exponentially large $`ae^{1/\rho J}`$, where $`a`$ is the lattice spacing and $`\rho `$ is the density of states at the Fermi level,E. S. Sørensen and I. Affleck (1996) with typical values on the order of 1 $`\mu `$m. Thus, it would increase as $`J`$ is decreased towards the physically relevant region $`J1`$, where we observe the peak at $`q_s=2k_F=\pi n`$. This might lead us to think that we should work with system sizes that are at least as large as the Kondo compensating cloud in order to observe features characteristic of a large FS. However, the fact that the peak at $`2k_F`$ becomes *sharper and more pronounced* as $`L`$ is increased (Fig. 3(b)) casts doubt on this naive expectation. Moreover, even if the true thermodynamic limit of the spin correlations do indeed oscillate at $`2k_F^{}=\pi \left(1+n\right)`$, our results show that perhaps *at the physically relevant distance scale* the important wave vector is actually $`2k_F=\pi n`$. For example, neutron scattering experiments are limited by the coherence length of the neutron beam and may not be able to probe large distances such as $`ae^{1/\rho J}`$. Other probes of the FS size, such as quantum oscillations, are limited by the electron mean free path, which would also have to exceed this length scale to access the asymptotic limit. Thus, our results put stringent limits on the observability of a large FS in heavy fermion systems. Besides, the presence of disorder and/or inelastic scattering may render the small FS size the only relevant length scale in real systems.

In conclusion, we have found clear signatures of a small Fermi surface in the spin correlation function of the one-dimensional Kondo lattice at small values of the Kondo coupling constant ($`J1`$). The discrepancy with previous results in the literature that had argued for a large Fermi surface in this system is ascribed to the presence of large edge perturbations introduced by the use of open boundary conditions. These perturbations are larger for large coupling constants ($`J1`$), which hinders the investigation of the Fermi surface size in this region. Even if the true asymptotic spin correlations are peaked at the large Fermi surface wave vector, the relevant oscillation period in many cases of interest may be set by the size of the small Fermi surface.

JCX is grateful to A. L. Malvezzi for providing some DMRG data for comparison. We also thank I. Affleck for constructive criticism and suggestions. This work was supported by the Brazilian agencies FAPESP (00/02802-7 and 01/00719-8) and CNPq (301222/97-5).


format: text

text:
# Spectral energy distributions of selfgravitating QSO discs

## 1 Introduction

The standard theoretical paradigm for the central engine in quasars and their radio-quiet kin (QSOs) is a viscous accretion disc surrounding a massive black hole. Direct evidence for accretion discs, such as double-peaked emission lines (Eracleous & Halpern, 1994), is important to seek but hard to come by. Perhaps the best reasons for belief in this paradigm are basic considerations of energy and angular momentum. Such discs are the most plausible astrophysical mechanism for converting rest mass to radiation with high efficiency, and if the relics of QSOs reside in galactic nuclei, then efficiencies $`10\%`$ are required (Soltan, 1982; Chokshi & Turner, 1992; Yu & Tremaine, 2002). On the other hand, the maximum specific angular momentum of a Kerr black hole, $`GM/c1.4M_8\mathrm{km}\mathrm{s}^1`$, is far less than that of most mass in galaxies. ($`M_8=M/10^8M_{}`$, where $`M`$ is the mass of the black hole.) Accreting gas must be separated from its angular momentum, and a viscous disc is the most natural mechanism.

Most of the binding energy of the disc is in its inner parts, so that the outer radius of a standard accretion disc is almost irrelevant to its bolometric luminosity. In contrast, most of the angular momentum is in the outer regions, so that the outer radius of the disc is sensitive to the initial angular momentum of the gas that feeds it. In a previous paper (Goodman, 2002, henceforth Paper I), it was argued that luminosity and angular momentum are coupled by the requirement that QSO discs be stable against their own self-gravity.

It is well known that gravitational stability is problematic in the outer parts of QSO discs (Shlosman & Begelman, 1987). The threat to the standard paradigm is that a strongly selfgravitating disc is likely to fragment completely into stars, leaving insufficient gas to fuel the QSO. As discussed in Paper I, viable solutions to this difficulty fall into several categories:

1. enhanced angular momentum transport, not necessarily by viscous processes but at rates corresponding to a viscosity parameter $`\alpha 1`$;
2. auxiliary heating in excess of what is provided by dissipating orbital energy, so as to reduce the gas density and selfgravity of the disc;
3. replacement of the outer disc with a very dense star cluster, whose collisional debris supply a disc of small radius and negligible selfgravity;
4. relatively low initial angular momentum for the gas, which therefore circularizes at small enough radius so as to avoid selfgravity.

A number of options in the first category were briefly considered in Paper I, including accretion driven by bars or global spiral waves (e.g. Shlosman & Begelman, 1989) or by magnetized winds (e.g. Blandford & Payne, 1982). But it was argued that none of these options is likely to achieve a supersonic accretion speed, and hence that discs cannot be stabilized much beyond one parsec by any of these mechanisms alone. The third category was tentatively rejected on the grounds that remnants of such star clusters are not observed in present day galactic nuclei. Recently, Pariev et al. (2002) have proposed discs whose thickness is supported primarily by magnetic pressure. If this is possible, it would imply lower disc densities on average, but the field might well squeeze the gas into dense clumps and thereby actually exacerbate selfgravity. At any rate, Pariev et al. (2002) do not apply their model to the cool outer regions beyond $`10^3R_\mathrm{S}`$.

The present paper will focus on some implications of the second category of solutions, which appear to be a natural compromise between a purely gaseous disc and a star cluster. It is reasonable to suppose that part of the gas fragments into stars, and that the energy released by nuclear fusion and other stellar processes (supernovae, stellar-mass black holes) may sufficiently heat the rest of the gas so as to prevent complete fragmentation (Collin & Zahn, 1999a, b). There is a great deal of evidence that such a feedback cycle operates in the discs of spiral galaxies on kiloparsec scales.

Observed SEDs of typical quasars differ markedly from classical theoretical predictions in which the disc is assumed to be geometrically thin, optically thick, steady, and heated solely by viscous dissipation (Pringle, 1981). To a first approximation, the typical SED is flat in a $`\lambda F_\lambda `$ plot over many decades in wavelength (Elvis et al., 1994). Relative to the classical predictions, there is excess emission at both X-ray and infrared wavelengths. The former is conventionally ascribed to comptonization in a hot corona at small radii (Shapiro et al., 1976), and the latter to passive reprocessing in warped or flared outer parts of the disc (Sanders et al., 1989).

We suppose that the infrared excess may be due to the energy inputs required to stabilize the outer disc against its own selfgravity. We obtain a lower bound on the auxiliary heating needed to stabilize the disc for given values of the macroscopic parameters: namely, the black hole mass ($`M`$), accretion rate ($`\dot{M}`$), and disc outer radius ($`r_{\mathrm{max}}`$). As shown in Paper I, and confirmed here with more realistic opacities, the inputs required for gravitational stability increase with the outer radius of the disc, $`r_{\mathrm{max}}`$. Paper I argued for an upper limit to $`r_{\mathrm{max}}`$ based on the energy available from plausible sources, such as fusion or accretion onto stellar-mass black holes. In this paper, we find limits to $`r_{\mathrm{max}}`$ from the SED.

It is not obvious that the auxiliary inputs should be completely thermalized, but if we assume this, then the SED can be predicted. The observed infrared emission of typical QSOs may be due largely to reprocessing of light emitted from the inner parts of the disc, as conventionally supposed. But by attributing all of the infrared light to the auxiliary energy inputs, and insisting that these inputs be sufficient for gravitational stability, we obtain bounds for $`r_{\mathrm{max}}`$. These bounds depend upon other parameters, especially the mass of the black hole, the accretion rate, and the viscosity parameter $`\alpha `$. We explore these dependencies. The meaning of $`r_{\mathrm{max}}`$ constrained by this method is the radius within which $`\dot{M}`$ is sensibly constant; obviously the disc can be extended indefinitely if the accretion rate and mass at large radii are sufficiently small. Hence our limits on $`r_{\mathrm{max}}`$ are best translated into upper limits on the initial angular momentum of the gas that is accreted.

The outline of our paper is as follows. §2 lays out the physical assumptions and governing equations for our disc models. These are the same as for the classical steady thin disc, with the one important exception that wherever the classical model would be gravitationally unstable, we invoke just enough auxiliary heating to stabilize it. A discussion of opacities becomes critically important, because the requirement of marginal stability fixes the density and temperature at the midplane (for given $`\dot{M}`$ and $`\alpha `$); the flux escaping from the disc, and hence the amount of auxiliary heating needed, then depend upon the optical depth. The computed radial structure of the disc is presented in §3 for parameters representative of bright quasars. Since our assumptions are no different from the classical ones at small radii, we emphasize the properties of our discs in the marginally selfgravitating region $`r10^3R_\mathrm{S}10^2\mathrm{pc}`$, and we compare our SEDs with those presented by Elvis et al. (1994) to obtain upper limits on the outer radius of the disc and the initial angular momentum of the gas. In the final section, we summarize our conclusions, issue the necessary caveats, and discuss directions for future research.

## 2 The selfgravitating accretion disc model

We start by considering the standard steady $`\alpha `$ disc. In this model, gas is accreted onto a central black hole of mass $`MM_8`$ in a steady, Keplerian, geometrically thin accretion disc at a constant rate $`\dot{M}`$. Accretion is driven by viscous mechanisms, but the exact mechanism (probably magnetorotationally-driven turbulence) is unimportant here. Newtonian equations will be used. These are not accurate near the inner edge of the disc, but the effects of selfgravity are at large radii where relativistic corrections are unimportant.

We assume that each annulus of the disc radiates as a blackbody with temperature $`T_{\mathrm{eff}}(r)`$. When accretion is the only source of thermal energy,

$$\sigma T_{\mathrm{eff}}^4=\frac{3}{8\pi }\dot{M}^{}\mathrm{\Omega }^2,$$
(1)
where $`\mathrm{\Omega }=(GM/r^3)^{1/2}`$, and we have defined for notational simplicity $`\dot{M}^{}=\dot{M}(1\sqrt{r_{\mathrm{min}}/r})`$ where $`r_{\mathrm{min}}`$ is the inner radius of the accretion disc. A zero-torque boundary condition has been assumed. In a relativistic thin-disc analysis, $`r_{\mathrm{min}}`$ would be the radius of the marginally stable orbit. In our newtonian framework, we set $`r_{\mathrm{min}}=R_\mathrm{S}/4ϵ`$, where $`R_\mathrm{S}=2GM/c^2`$ is the Schwarzschild radius and $`ϵ0.1`$ is the assumed radiative efficiency of accretion.

An $`\alpha `$ prescription for the viscosity is adopted, $`\nu =\alpha c_\mathrm{s}h\beta ^b`$, where $`c_\mathrm{s}=\sqrt{p/\rho }`$ is the isothermal sound speed at the disc midplane, $`h=c_\mathrm{s}/\mathrm{\Omega }`$ the half thickness, and $`\beta =p_{\mathrm{gas}}/p`$ is the ratio of gas pressure to total pressure at the midplane. The viscosity is related to the accretion rate by

$$\beta ^bc_\mathrm{s}^2\mathrm{\Sigma }=\frac{\dot{M}^{}\mathrm{\Omega }}{3\pi \alpha },$$
(2)
where $`\mathrm{\Sigma }=_{\mathrm{}}^+\mathrm{}\rho 𝑑z`$ is the surface density at a given radius. The parameter $`b`$ is a switch that can be either $`1`$, so that viscosity is proportional to gas pressure, or $`0`$, so that viscosity is proportional to total pressure. Magnetic pressure is assumed to contribute negligibly to the thickness of the disc, even if magnetic stresses dominate angular momentum transport.

Because it will turn out that the optically thick assumption fails at large radii, we now turn to a study of the relationship between the midplane temperature $`T`$ and $`T_{\mathrm{eff}}`$. We assume radiative energy transport in the vertical ($`z`$) direction and write $`\tau =_0^{\mathrm{}}\kappa \rho 𝑑z`$ for the optical depth at the midplane. In the optically thick limit $`\tau 1`$, the diffusive approximation is applicable so that $`T^4\tau T_{\mathrm{eff}}^4`$, with a numerical factor of order unity that depends upon the vertical variation of opacity $`\kappa `$ and the vertical distribution of viscous dissipation. The latter is not yet predictable from theory, despite major advances in understanding magnetohydrodynamic disc turbulence. If one assumes constant dissipation per unit optical depth, then

$$T^4\left(\frac{3}{8}\tau +\frac{1}{2}\right)T_{\mathrm{eff}}^4,$$

where the “$`+1/2`$” follows from the Eddington approximation to the boundary condition at the photosphere. In the opposite limit of low optical depth, the radiative flux from one side of the disc scales as $`\sigma T_{\mathrm{eff}}^4\tau \sigma T^4`$. Again a numerical factor occurs. One should distinguish between absorption and scattering opacity, but at the large accretion rates of interest to us, optically thin conditions tend to occur where the opacity is primarily absorptive. For a vertically isothermal disc, $`\sigma T_{\mathrm{eff}}^4=4\tau \sigma T^4`$ when $`\tau 1`$. Simple interpolation between the optically thin and optically thick limits yields eq. (4) below, which satisfies the physical requirement that $`T_{\mathrm{eff}}T`$ for any $`\tau `$ (in fact, the minimum of $`T^4/T_{\mathrm{eff}}^4`$ is $`1.112`$ and occurs when $`\tau =\sqrt{2/3}`$).

We represent the effective radiation pressure by eq. (7). In the optically thick case, in view of eq. (4), this reduces to the usual $`p_{\mathrm{rad}}=aT^4/3`$ relation, but in the optically thin case, $`p_{\mathrm{rad}}2\tau ^2\sigma T^4/c`$: one power of $`\tau `$ represents the inefficiency with which photons are radiated, and the second is the fraction of these photons that transfer their momentum to the gas, thereby helping to support the thickness of the disc. The numerical factor is correct for an isothermal slab of midplane optical depth $`\tau 1`$.

We also must know the opacity in the disc (equation 13). We use the opacity tables of Iglesias & Rogers (1996) for high temperatures and those of Alexander & Ferguson (1994) for low temperatures, with $`X=.70`$, $`Z=.03`$. Figure 1 illustrates the opacity as a function of density and temperature. The two opacity tables overlap in the region $`3.75<\mathrm{log}T<4.1`$, and Figure 1 demonstrates their agreement. Additionally, the loci of two disc models are projected onto the plot. Extrapolation of the opacity tables in the low density regime is necessary. This is done by simply assigning $`\kappa `$ its value at the low-density edge of the opacity table (using the same value of $`T`$). This extrapolation is presumably correct for high temperatures ($`T10^4\mathrm{K}`$) where the opacity is dominated by electron scattering. At lower temperatures, the extrapolation is less certain, but we feel it is reasonable as a first approximation. Also, for $`\mathrm{log}T<3`$, $`\kappa `$ is assigned its value at the low-temperature edge of the opacity table, which is approximately $`10^{0.76}\mathrm{cm}^2\mathrm{g}^1`$.

In summary, our governing equations for the disc are as follows (we take $`m=.62m_H`$ as the mean molecular mass):

$`\sigma T_{\mathrm{eff}}^4`$ $`=`$ $`{\displaystyle \frac{3}{8\pi }}\dot{M}^{}\mathrm{\Omega }^2`$ (3)
$`T^4`$ $`=`$ $`\left({\displaystyle \frac{3}{8}}\tau +{\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{4\tau }}\right)T_{\mathrm{eff}}^4`$ (4)
$`\tau `$ $`=`$ $`{\displaystyle \frac{\kappa \mathrm{\Sigma }}{2}}`$ (5)
$`\beta ^bc_\mathrm{s}^2\mathrm{\Sigma }`$ $`=`$ $`{\displaystyle \frac{\dot{M}^{}\mathrm{\Omega }}{3\pi \alpha }}`$ (6)
$`p_{\mathrm{rad}}`$ $`=`$ $`{\displaystyle \frac{\tau \sigma }{2c}}T_{\mathrm{eff}}^4`$ (7)
$`p_{\mathrm{gas}}`$ $`=`$ $`{\displaystyle \frac{\rho kT}{m}}`$ (8)
$`\beta `$ $`=`$ $`{\displaystyle \frac{p_{\mathrm{gas}}}{p_{\mathrm{gas}}+p_{\mathrm{rad}}}}`$ (9)
$`\mathrm{\Sigma }`$ $`=`$ $`2\rho h`$ (10)
$`h`$ $`=`$ $`{\displaystyle \frac{c_\mathrm{s}}{\mathrm{\Omega }}}`$ (11)
$`c_\mathrm{s}^2`$ $`=`$ $`{\displaystyle \frac{p_{\mathrm{gas}}+p_{\mathrm{rad}}}{\rho }}`$ (12)
$`\kappa `$ $`=`$ $`\kappa (\rho ,T)`$ (13)
For an accretion disc with given parameters $`M`$, $`\dot{M}`$, $`\alpha `$, and $`b`$, we can solve these eleven equations for the eleven unknowns: $`T_{\mathrm{eff}}`$, $`T`$, $`\tau `$, $`\mathrm{\Sigma }`$, $`\beta `$, $`c_\mathrm{s}`$, $`p_{\mathrm{rad}}`$, $`p_{\mathrm{gas}}`$, $`\rho `$, $`h`$, and $`\kappa `$ at every radius $`r`$, which is related simply to $`\mathrm{\Omega }`$.

The outer parts of the accretion disc will be prone to selfgravity if Toomre’s stability parameter,

$$Q\frac{c_\mathrm{s}\mathrm{\Omega }}{\pi G\mathrm{\Sigma }}\frac{\mathrm{\Omega }^2}{2\pi G\rho }$$
(14)
is $`1`$. In our model, we assume there is some feedback mechanism, perhaps star formation, that supplies just enough additional heat to the disc as to prevent $`Q`$ falling below a minimum value $`Q_{\mathrm{min}}1`$. The energy source for these auxiliary inputs must be something other than the orbital energy of the gas, for example, nuclear fusion in stars. The quantity of these inputs is fixed by the properties of the disc: in particular the opacities, since the more easily radiation escapes, the more the disc must be heated to maintain stability. The nature of the auxiliary sources is not important to the spectral energy distribution as long as their inputs are completely thermalized.

When $`Q=Q_{\mathrm{min}}`$, eq. (3), which assumes the disc to be heated by release of orbital energy only, is no longer appropriate; we replace it with

$$\rho =\frac{\mathrm{\Omega }^2}{2\pi GQ_{\mathrm{min}}}$$
(15)
and $`T_{\mathrm{eff}}`$ must be obtained from equation 4.

Once $`T_{\mathrm{eff}}`$ is known as a function of $`r`$, we compute a theoretical spectral energy distribution for the disc as a weighted sum of Planck functions:

$$L_\lambda =\frac{2hc^2}{\lambda ^5}_{r_{\mathrm{min}}}^{r_{\mathrm{max}}}2\pi r𝑑r\frac{1}{\mathrm{exp}[hc/\lambda kT_{\mathrm{eff}}(r)]1}$$
(16)
This formula does not allow for differences between effective and color temperature. Such differences are most likely in those parts of the disc whose emission is dominated by passive reprocessing (Sanders et al., 1989) or are optically thin (to absorption). In fact we find that much of the auxiliary heat input occurs in regions that are marginally optically thick, since these regions cool most efficiently. We feel that eq. (16) is adequate for a first investigation: it may somewhat overestimate the wavelength of the infrared peak associated with the auxiliary heating, but the area under the peak should be roughly correct.

## 3 Results

We define $`l_\mathrm{E}L_0/L_\mathrm{E}`$, where $`L_0ϵ\dot{M}c^2`$ is the luminosity corresponding to a nonselfgravitating (i.e., with no auxiliary heating) disc with accretion rate $`\dot{M}`$, and $`L_\mathrm{E}`$ is the Eddington luminosity. In our models we always set $`ϵ=.1`$ and $`Q_{\mathrm{min}}=1`$.

The radial structure of a typical accretion disc ($`M_8=1`$, $`l_\mathrm{E}=.5`$, $`\alpha =.01`$) is presented in Figure 2 for $`b=0`$. The outer radius of $`r_{\mathrm{max}}=10^7R_\mathrm{S}`$ will be shown to be unrealistic, but this model illustrates several key regimes in the disc.

The inner region ($`2.5r/R_\mathrm{S}10^2`$) is radiation dominated ($`\beta 1`$) and the temperature is high enough to ensure that the opacity is dominated by electron scattering, hence constant. In view of equation 1, $`T_{\mathrm{eff}}r^{3/4}`$. When $`Q`$ would otherwise drop below $`1`$ (here, at $`r10^3R_\mathrm{S}`$), we hold $`Q`$ constant, which ensures that $`\rho r^3`$. In this auxiliary heating region, if $`\kappa =\mathrm{constant}`$, then $`T_{\mathrm{eff}}r^{3/8}`$. After the disc becomes optically thin ($`\tau <1`$) we find multiple solutions. The one with the least optical depth corresponds to the continuous line in all panels; it reaches constant $`T`$ and $`\kappa `$. There is an optically thick solution with $`Tr^{3/4}`$. Since the region in which we find multiple solutions is so dependent on the extrapolation of opacity into low-density and low-temperature regimes, the interpretation must be approached with caution. However, the presence of multiple solutions raises the possibility of thermal instability and time-dependence (see, for example, Siemiginowska et al., 1996). Note that in our solutions, $`\dot{M}`$ is constrained to be constant throughout the disc. Several of the panels in Fig. 2 contain only one curve in the region of multiple solutions ($`r10^5\mathrm{pc}`$): this is because $`\rho `$ and $`c_\mathrm{s}`$ are uniquely determined by $`r`$, $`\dot{M}`$, $`\alpha `$, and $`Q`$ *via* eqs. (15) and (6) when $`b=0`$.

For $`b=1`$, the effective temperature has the same $`T_{\mathrm{eff}}r^{3/4}`$ behavior in the inner part of the disc, but differs from the $`b=0`$ case in the auxiliary heating region. For $`\kappa \mathrm{constant}`$ in this region, in the optically thick solution, $`Tr^{1/2}\&T_{\mathrm{eff}}r^{1/4}`$, and in the optically thin case, $`T\text{const.}\&T_{\mathrm{eff}}r^{3/8}`$.

The calculated SED of this typical model is presented in Figure 3, for outer radii of $`r_{\mathrm{max}}=10^3`$, $`10^4`$, $`10^5`$, $`10^6`$, and $`10^7R_\mathrm{S}`$. (All of our SEDs assume that the disc is viewed face-on.) Since $`T_{\mathrm{eff}}`$ may have multiple solutions in certain regions of the disc, the SEDs are plotted as bounded regions with the top (bottom) boundary corresponding to the largest (least) $`T_{\mathrm{eff}}`$. Also plotted is the arbitrarily normalized mean energy distribution (MED) from 47 observed quasars (Elvis et al., 1994). The authors caution that the dispersion about this mean is large, of order “one decade for both the infrared and ultraviolet components when the MED is normalized at \[$`\lambda =1.25\mu \mathrm{m}`$\].” The $`r_{\mathrm{max}}=10^3R_\mathrm{S}`$ case has an expected $`\lambda f_\lambda \lambda ^{4/3}`$ behavior since selfgravity is not an issue. Relativistic effects (Laor & Netzer, 1989; Sun & Malkan, 1989) only influence the SED at and blueward of the short wavelength peak, since this emission comes from the innermost parts of the disc.

Figure 4 shows the dependence of the infrared bump on the mass of the central black hole. (Note that the wavelength corresponding to the blue bump has an expected $`M^{1/4}`$ dependence.) The Eddington ratio as defined above is kept constant at $`l_\mathrm{E}=.5`$: note that this parameter defines $`\dot{M}`$ rather than the total luminosity of the disc, which includes the auxiliary heat inputs. In the top two panels, $`r_{\mathrm{max}}`$ is set to $`10^5R_\mathrm{S}`$, and the relative prominence of the infrared bump grows with $`M`$ because the reduced radius at which $`Q_{\mathrm{min}}=1`$ decreases with $`M`$ (Paper I). The bottom two panels have $`r_{\mathrm{max}}`$ fixed at $`1\mathrm{pc}`$; here larger values of $`M`$ tend to reduce the infrared bump. This dependence is illustrated in Figure 7. Solutions for the disc are considered incompatible with the observations if $`\lambda f_\lambda `$ at the maximum of the infrared bump is more that ten times that of the blue bump.

Figure 5(a) shows that while the blue bump increases with mass accretion rate, the energy output of the infrared bump is hardly affected, for the $`b=0`$ case. This invariance is caused by two competing factors as $`\dot{M}`$ is increased: $`T_{\mathrm{eff}}`$ is higher in the selfgravitating region, but because $`T`$ is also higher, the onset of the optically thin regime, and thus a sharp increase in $`T_{\mathrm{eff}}`$ as shown in Figure 2, is delayed. Similar considerations explain the behavior of the infrared bump in Figure 5(b), for $`b=1`$. For low accretion rates, the midplane temperature drops low enough that the disc becomes optically thin before the end of the disc (at $`r_{\mathrm{max}}=10^5R_\mathrm{S}`$), which drives up $`T_{\mathrm{eff}}`$.

In Figure 6, the same argument explains the behavior of the infrared bump in both the $`b=0`$ and $`b=1`$ cases. The onset of the optically thin selfgravitating region is at smaller radii with increasing $`\alpha `$ for both cases, but $`T_{\mathrm{eff}}`$ is higher in the optically thick selfgravitating region for $`b=1`$.

It is clear that additional heating in the outer parts of our accretion discs results in a second bump in the SED. Assuming that this outer-disc contribution to the SED can be no more prominent than infrared emission in the observed MED, we can already place tight constraints on the outer radius of these disc models: $`r_{\mathrm{max}}10^410^5R_\mathrm{S}0.11\mathrm{pc}`$.

## 4 Summary and Discussion

We have estimated spectral energy distributions (SEDs) of bright QSOs using standard assumptions, with one addition: where the disc would otherwise be gravitationally unstable, we have postulated additional sources of heat, other than release of orbital energy by accretion, just sufficient to maintain gravitational stability. These sources become necessary beyond $`10^3R_\mathrm{S}`$ for typical parameters. Assuming their energy inputs to be locally and completely thermalized, we have calculated their contribution to the SED and luminosity of the disc, which occurs primarily in the red and near infrared. The larger the disc, the more auxiliary heating is required. For typical black-hole masses and accretion rates inferred from the blue bump, the auxiliary inputs actually exceed the power derived from accretion if the disc extends beyond $`10^410^5R_\mathrm{S}`$, or about one parsec. This would be incompatible with the typical SED of bright QSOs, which is approximately flat in $`\lambda F_\lambda `$.

Paper I placed similar limits on $`r_{\mathrm{max}}`$ from energetic arguments. It was assumed that the stars or small black holes that heat the disc also form within the disc, and that the mass in these objects is at most comparable to that of the disc. The present limits do not rely on these assumptions. That is to say, if the disc were heated by a much larger mass in stars, the limits of Paper I could be relaxed, but those based on the SED would still apply.

It has been assumed that the accretion rate is constant with radius, so that $`\dot{M}`$ at large radii can be derived from the luminosity in the blue bump (which has no contribution from the auxiliary sources). If $`\dot{M}`$ at $`r1\mathrm{pc}`$ is several orders of magnitude less than it is at $`r10^3R_\mathrm{S}`$, then the disc could be much more extensive than we have supposed. The Stefan-Boltzmann law implies that the observed luminosity at $`10\mu `$ must in any case come from $`r2(\lambda L_\lambda /10^{46}\text{ erg s}^1)^{1/2}\mathrm{pc}`$, but this does not require that the surface density at that distance is as high as in a constant-$`\dot{M}`$ disc. Therefore, it is useful to rephrase the limits in terms of the initial angular momentum of the gas supplied to the disc. At least on a time average, $`\dot{M}`$ should be constant inside the radius corresponding to the initial angular momentum. The relationship is

$$J_0660\left(r_{\mathrm{pc}}M_8\right)^{1/2}\mathrm{pc}\mathrm{km}\mathrm{s}^1.$$
(17)
This is quite small compared with the product of virial velocity ($`300\mathrm{km}\mathrm{s}^1`$) and scale size ($`1\mathrm{kpc}`$) of QSO hosts, and it may be important to ask where gas with such low angular momentum comes from. Furthermore, the mass of a gravitationally stable disc that obeys our constraints on the SED is generally much less than that of the black hole, so that the disc must be replenished many times over to grow the black hole by accretion.

Our treatment of the disc is highly simplified and certainly crude compared with many past efforts. We have ignored relativistic effects, adopted a one-zone model for the vertical structure, used Rosseland mean opacities without distinguishing between scattering and absorption, taken a constant molecular weight, and represented angular-momentum transport by the usual viscous $`\alpha `$ prescription. We feel that these simplifications are justified by the strong dependence of the auxiliary inputs and the SED on the outer radius of the disc. A more detailed treatment of the physics that one actually understands seems unlikely to change our conclusions concerning $`r_{\mathrm{max}}`$. Radical enhancements in transport (equivalent to $`\alpha 1`$), or magnetic pressures $`p_{\mathrm{gas}}`$ (as suggested by Pariev et al., 2002) could make some difference, but these are not yet understood.

Despite what has just been said, a more elaborate treatment of vertical structure and radiative transfer might point to a redistribution of the auxiliary energy inputs in wavelength, if not their total contribution to the disc luminosity. Unfortunately, one has no reliable predictions for the vertical distribution of purely viscous heating, much less of the auxiliary sources postulated here, which will limit the credibility of detailed vertical models.

Even if the actively accreting parts of QSO discs are smaller than $`0.1\mathrm{pc}`$, it is still possible that selfgravity is important in them, and therefore that they form stars. Quiescent galactic nuclei with black holes, even in early-type galaxies, often show kinematic evidence of compact stellar discs (Gebhardt et al., 2000; Bower et al., 2001; de Zeeuw et al., 2002). Nuclear starbursts (albeit on scales $`10^2\mathrm{pc}`$) are usually accompanied by AGN activity (Sanders, 1999; Heckman, 1999). The black hole in our own Galaxy, though not an AGN and estimated to have a very low accretion rate (Quataert et al., 1999), is surrounded by what appear to be young high-mass stars at $`r0.1\mathrm{pc}`$ (Krabbe et al., 1995). More theoretical attention should be paid to star formation in these extreme environments, where the densities, temperatures, and tidal fields are much higher than in normal giant molecular clouds.

We thank Iskra Strateva and Jonathan Tan for helpful discussions.


text:
# QCD phase diagram for small densities from simulations at imaginary 𝜇

## 1 INTRODUCTION

In view of heavy ion collision experiments a pressing task for lattice QCD is to map out the deconfinement transition in the $`(\mu ,T)`$-plane. The complexity of the fermion determinant at finite baryon density renders standard Monte Carlo techniques impossible . One way of circumventing this ‘sign problem’ are multi-dimensional reweighting techniques, and a phase diagram has been presented last year . Reweighting generally breaks down for large volumes and/or densities, and one difficulty of the approach is to know when that happens. In  both reweighting factor and observables were expanded in a Taylor series in $`\mu /T`$, and the first coefficients were measured . This method allows to simulate larger volumes, but a priori nothing is known about the error introduced by truncating the series.

Here we present an alternative method avoiding reweighting altogether. We simulate at imaginary $`\mu `$ for which there is no sign problem and arbitrarily large volumes are feasible. Observables are fitted by truncated Taylor series in $`\mu /T`$, keeping full control over the associated systematic error. The series may then be continued to real $`\mu `$, provided the observable is analytic . This strategy has been successfully tested for screening masses , here we extend it to the critical line itself. A detailed account of this work is given in .

## 2 QCD AT IMAGINARY $`\mu `$

Finite density QCD is symmetric under reflection of $`\mu `$, $`Z(\mu )=Z(\mu )`$, for both real and imaginary $`\mu `$. For imaginary $`\mu =\mathrm{i}\mu _I`$ it is moreover periodic, $`Z(\mu _I)=Z(\mu _I+2\pi /3)`$ . Hence, certain shifts in $`\mu _I`$ are equivalent to Z(3) transformations to different vacua. Consequently, Z(3) transitions occur at the critical values $`(\mu /T)^c=(k+1/2)2\pi /3`$. The Z(3) sectors can be identified by monitoring the phase of the Polyakov loop. Ref.  predicts these transitions to be of first order for $`T>T_c`$ and crossover for $`T<T_c`$, with the critical deconfinement temperature $`T_c`$. Schematically, the phase diagram is shown in Fig. 1, which is periodically repeated for larger $`\mu _I`$.

Our simulations support this picture. Fig. 2 shows the phase of the Polyakov loop at $`(\mu /T)^c`$ for various lattice couplings/temperatures. A smooth distribution is found for low temperatures, whereas at high temperatures the two-state distribution signals a first order phase transition. Similar findings are reported for four flavors in . The existence of a Z(3) transition with its endpoint joining on the deconfinement transition implies that the critical line of the latter becomes non-analytic at that point, and hence only $`\mu _I\pi /3`$ can be continued to real $`\mu `$.

Note that the deconfinement line as well as all observables even under reflections $`\mu \mu `$ are symmetric around the Z(3) transition. This is a feature that is not reproduced by reweighting from a Monte Carlo ensemble generated at $`\mu =0`$, since such an ensemble does not sufficiently probe the other Z(3) sectors .

## 3 CONTINUATION OF THE DECONFINEMENT LINE

One method to locate a phase transition is to look for the peak $`\chi _{max}=\chi (\mu _c,\beta _c)`$ of the susceptibilities $`\chi =VN_t(𝒪𝒪)^2`$, where we have used the plaquette, the chiral condensate and the Polyakov loop as operators $`𝒪(x)`$. Since simulations are always performed on finite volumes, susceptibilities are analytic functions over the whole parameter space of the theory. Non-analyticities associated with phase transitions develop only in the thermodynamic limit.

For a given finite volume, the pseudo-critical line is then found by locating $`\chi _{max}`$ for every $`\mu `$,

$$\frac{\chi }{\beta }|_{\beta _c}=0\frac{^2\chi }{\beta ^2}|_{\beta _c}<0.$$
(1)
This is an implicit definition of the critical coupling $`\beta _c(\mu )`$. For analytic $`\chi (\mu ,\beta )`$ the implicit function theorem implies that also $`\beta _c(\mu )`$ is analytic for all $`\mu `$. Moreover, the symmetry considerations of the last section ensure $`\beta _c(\mu )=\beta _c(\mu )`$, so that we can expand in $`\mu /T=aN_t\mu `$ to obtain

$$\beta _c(\mu )=\underset{n}{}c_n(a\mu )^{2n}.$$
(2)
As long as, for imaginary $`\mu =\mathrm{i}\mu _I`$, the non-perturbative data are well described by a few coefficients in this series, analytic continuation to real $`\mu `$ simply proceeds by the operation $`\mu _I\mathrm{i}\mu _I`$. This provides us with the location of the pseudo-critical line on a finite volume.

If there is a first order deconfinement transition ending in a critical point, then, as the volume is increased, $`\beta _c(\mu )`$ approaches the line of phase transitions, while remaining pseudo-critical in a crossover regime. In principle, the nature of the line can be determined from the way the thermodynamic limit is approached, $`(\beta _c(V)\beta _c(\mathrm{}))V^\sigma `$, with $`\sigma =1`$ or $`0<\sigma <1`$ for a first or second order transition, respectively, and $`\sigma =0`$ for a crossover. Our first simulations are on a single volume and this question is postponed to future work.

Fig.3 shows a sequence of susceptibilities for the modulus of the Polyakov loop, with similar pictures for the plaquette and the chiral condensate. All three observables give statistically compatible values $`\beta _c(\mu _I)`$. Note that the peak is moving to larger $`\beta _c`$ for growing $`\mu _I`$, as one would expect for imaginary $`\mu `$, cf. Eq. 2 and Fig. 1.

Fig. 4 shows $`\beta _c`$ for different values of $`\mu _I`$. The data perfectly reproduce the expected symmetry of the partition function about the first Z(3) transition point. The lines in the figure show fits of a quadratic and a quartic polynomial to the data, with $`\chi ^2/\mathrm{dof}`$0.6 and 0.7, respectively. The coefficient of the quartic term is found to be zero within errors, and hence the two fits are statistically compatible: up to $`\mu _I/T=\pi /3`$, and for the quark flavors and masses under consideration, the critical line is well described by the quadratic term only.

After analytic continuation the result is transformed into physical units by means of the perturbative two-loop expression for the scale parameter $`a\mathrm{\Lambda }_L`$ and using $`T_c(\mu =0)`$ from the literature  to set the physical scale. We thus arrive at the phase diagram of Fig. 5, with the pseudo-critical line given by ($`\mu _B=3\mu `$)

$$\frac{T_c(\mu _B)}{T_c(\mu _B=0)}=10.00563(38)\left(\frac{\mu _B}{T}\right)^2.$$
(3)
This result is consistent with those obtained by the reweighting approaches , for which it provides a crucial, controlled check. An important task ahead is to perform a finite size scaling analysis of the continued critical line, in order to determine the order of the transition as well as to locate a possible critical point.


text:
# Separating cosmic shear from intrinsic galaxy alignments: correlation function tomography

## 1 Introduction

The tidal gravitational field of mass inhomogeneities distorts the images of distant galaxies, resulting in correlations in their observed ellipticities. This cosmological weak lensing signal, or cosmic shear, depends upon cosmological parameters and the matter power spectrum (Blandford et al. 1991; Miralda-Escudé 1991; Kaiser 1992). In 2000, four teams announced the first detections of cosmic shear (Bacon et al. 2000; Kaiser et al. 2000; van Waerbeke et al. 2000; Wittman et al. 2000; Maoli et al. 2001), and more recently measurements at arcminute scales have been made using the HST (Hämmerle et al. 2002; Refregier et al. 2002). Interesting constraints have already been placed on the matter power spectrum normalisation $`\sigma _8`$, and cosmic shear is also particularly sensitive to the matter density parameter $`\mathrm{\Omega }_\mathrm{m}`$, the power spectrum shape parameter $`\mathrm{\Gamma }`$ and the source redshift distribution (e.g. van Waerbeke et al. 2002a; Hoekstra et al. 2002). Future multi-colour surveys will cover hundreds of square degrees, and have the potential to place tight constraints on cosmological parameters particularly when combined with results from the CMB, SNIa and galaxy surveys (Mellier et al. 2002; van Waerbeke et al. 2002b). For example, van Waerbeke et al. (2002a) compared their constraints on $`\mathrm{\Omega }_\mathrm{m}`$ and $`\sigma _8`$ from lensing with those of Lahav et al. (2002) from the CMB, noting their near orthogonality.

Various statistical measures of the cosmic shear have been suggested; here we focus on the two-point shear correlation function $`\xi _+`$ (hereafter denoted by $`\xi `$), which is convenient since it is insensitive to gaps in the data field, unlike integrated measures such as the aperture mass statistic $`_{\mathrm{ap}}`$ (e.g. Schneider et al. 1998).

A possible systematic contaminant of the lensing correlation function $`\xi ^\mathrm{L}`$ is intrinsic alignment, which may arise during the galaxy formation process. This has been subject to numerical, analytic and observational studies \[e.g. Croft & Metzler 2000; Heavens et al. 2000 (HRH); Crittenden et al. 2001; Catelan et al. 2001; Mackey et al. 2002; Brown et al. 2002; Jing 2002; Hui & Zhang 2002\], where amplitude estimates span a few orders of magnitude due to differences in the mechanism assumed to be responsible, and the type of galaxy considered. Nevertheless, these studies agree that the intrinsic correlation signal $`\xi ^\mathrm{I}`$ can dominate the lensing signal for surveys with $`z\begin{array}{c}<\hfill \\ \hfill \end{array}0.5`$.

Any correlation in ellipticities due to intrinsic alignment only arises from physically close galaxy pairs, whereas $`\xi ^\mathrm{L}`$ is sensitive to the integrated effect of the density fluctuations out to the redshift of the nearer galaxy. It has been shown that photometric redshift information could be used to suppress $`\xi ^\mathrm{I}`$, by downweighting or ignoring galaxy pairs at approximately the same redshift (King & Schneider 2002), or by downweighting nearby pairs and subtracting a model of the intrinsic alignment signal from the observed ellipticity correlation function (Heymans & Heavens 2002).

Motivated by the fact that intrinsic galaxy alignment is not yet well understood, we present a new method to isolate the intrinsic and lensing-induced components of the galaxy ellipticity correlation function. This method assumes that photometric redshift information is available, so that the correlation function can be measured between different redshift slices. However, no specific model for intrinsic ellipticity correlation (for instance its correlation length or redshift evolution) needs to be adopted. In the next section we outline the method and in Sect. 3 we present some results in the context of a possible future survey. We discuss the results in Sect. 4.

## 2 General method and specific assumptions

In this section we outline a method to separate and extract the intrinsic and lensing components of the galaxy ellipticity correlation function. The general method is described in Sect. 2.1, and in Sect. 2.2 and Sect. 2.3 we state the assumptions particular to this work.

### 2.1 Method

The ellipticity correlation function for galaxies with angular separation $`\theta `$ and at true redshifts $`z_i`$, $`z_j`$, is composed of a lensing and an intrinsic signal

$$\xi (\theta ,z_i,z_j)=\xi ^\mathrm{L}(\theta ,z_i,z_j)+\xi ^\mathrm{I}(\theta ,z_i,z_j).$$
(1)
As noted in Sect. 1, the origin and behaviour of $`\xi ^\mathrm{I}`$ is not yet well understood. $`\xi ^\mathrm{L}(\theta ,z_i,z_j)`$ is related to the 3-dimensional matter power spectrum $`P_\delta `$ through

$`\xi ^\mathrm{L}(\theta ,z_i,z_j)={\displaystyle \frac{9H_0^4\mathrm{\Omega }_\mathrm{m}^2}{4c^4}}{\displaystyle _0^{\mathrm{min}[w_i,w_j]}}{\displaystyle \frac{\mathrm{d}w}{a^2(w)}}`$
$`\times R(w,w_i)R(w,w_j){\displaystyle \frac{\mathrm{d}\mathrm{}\mathrm{}}{(2\pi )}P_\delta (\frac{\mathrm{}}{f(w)},w)\mathrm{J}_0(\mathrm{}\theta )},`$ (2)
where $`H_0`$ and $`\mathrm{\Omega }_\mathrm{m}`$ are the values of the Hubble parameter and matter density parameter at the present epoch, and $`a(w)`$ is the scale factor at comoving distance $`w`$, normalised such that $`a(0)=1`$ today. J0 is the $`0`$-th order Bessel function of the first kind, and $`\mathrm{}`$ is the angular wave-vector. We denote the comoving distance at $`z_i`$ by $`w_i`$, and the function $`f(w)`$ is the comoving angular diameter distance, which depends on the spatial curvature $`K`$:

$$f(w)=\left\{\begin{array}{cc}K^{1/2}\mathrm{sin}(K^{1/2}w)\hfill & (K>0)\hfill \\ w\hfill & (K=0)\hfill \\ (K)^{1/2}\mathrm{sinh}[(K)^{1/2}w]\hfill & (K<0)\hfill \end{array}\right\}.$$
(3)
The function $`R(w,w^{})=f(w^{}w)/f(w^{})`$ is the ratio of the angular diameter distance of a source at comoving distance $`w^{}`$ seen from a distance $`w`$, to that seen from $`w=0`$.

Next, we account for the availability of photometric redshift estimates rather than spectroscopic ones. The galaxy ellipticity correlation function becomes

$$\overline{\xi }(\theta ,\overline{z}_i,\overline{z}_j)=dz_idz_jp(z_i,z_j|\overline{z}_i,\overline{z}_j,\theta )\xi (\theta ,z_i,z_j),$$
(4)
where $`p(z_i,z_j|\overline{z}_i,\overline{z}_j,\theta )`$ is the probability to have true redshifts $`z_i`$ and $`z_j`$ given photometric estimates $`\overline{z}_i`$ and $`\overline{z}_j`$. This is given by

$`p(z_i,z_j|\overline{z}_i,\overline{z}_j,\theta )=`$ (5)
$`{\displaystyle \frac{p(z_i|\overline{z}_i)p(z_j|\overline{z}_j)[1+\xi _{\mathrm{gg}}(r)]}{dz_idz_jp(z_i|\overline{z}_i)p(z_j|\overline{z}_j)[1+\xi _{\mathrm{gg}}(r)]}},`$
where $`\xi _{\mathrm{gg}}`$ is the galaxy spatial correlation function, which may also include redshift dependent evolution, and $`r`$ is the comoving separation of the galaxies in a pair.

We now assume that $`\overline{\xi }(\theta ,\overline{z}_i,\overline{z}_j)`$ is available on a 3-dimensional grid of $`N_K`$ angular separation bins of width $`\mathrm{\Delta }\theta `$ centred on $`\theta _K`$ (index $`K`$), and $`N_Z`$ photometric redshift bins of width $`\mathrm{\Delta }z`$ centred on each of $`\overline{z}_i`$ (index $`I`$) and $`\overline{z}_j`$ (index $`J`$). This could either correspond to an observed signal $`\overline{\xi }_{IJK}^{\mathrm{obs}}`$, or to a theoretical prediction $`\overline{\xi }^{\mathrm{mod}}_{IJK}`$ which we want to compare with the observed signal.

We assume that both the lensing and intrinsic correlations can be written in terms of sets of template functions $`A_n`$ and $`B_n`$

$`\overline{\xi ^\mathrm{L}}(\theta ,\overline{z}_i,\overline{z}_j)={\displaystyle \underset{n=1}{\overset{\mathrm{N}_\mathrm{L}}{}}}a_nA_n(\theta ,\overline{z}_i,\overline{z}_j),`$ (6)
$`\overline{\xi ^\mathrm{I}}(\theta ,\overline{z}_i,\overline{z}_j)={\displaystyle \underset{n=1}{\overset{\mathrm{N}_\mathrm{I}}{}}}b_nB_n(\theta ,\overline{z}_i,\overline{z}_j),`$
where $`a_n`$ and $`b_n`$ are the amplitudes of the $`n`$-th lensing and $`n`$-th intrinsic template functions. The template functions are fairly arbitrary, and extra functions can be added as required, to span the range of plausible models. We describe our choice of models in Sect. 2.3 below.

A suitably chosen single index $`m`$ identifies correlations between bins with redshift indices $`I`$,$`J`$ and angular separation index $`K`$. In total, there are $`N_M=N_Z(N_Z+1)N_K/2`$ such independent measurements. The total of $`N=N_L+N_I`$ gridded template models for the correlation functions ($`A_1\mathrm{}A_{N_L},B_1\mathrm{}B_{N_I}`$) can be written as an $`N_M\times N`$ so-called design matrix $``$, and their amplitudes ($`a_1\mathrm{}a_{N_L},b_1\mathrm{}b_{N_I}`$) as an $`N`$-dimensional column vector $`𝒢`$ so that

$$\overline{\xi }^{\mathrm{mod}}_m=_{mn}𝒢_n.$$
(7)
Our aim is to recover $`\overline{\xi }_m^{\mathrm{obs}}`$ in terms of the template functions. Using the method of least squares, the best-fit estimates of the $`𝒢_n`$ are those values $`\widehat{𝒢}_n`$ which minimise

$$S=(\overline{\xi }^{\mathrm{obs}}𝒢)𝒞^1(\overline{\xi }^{\mathrm{obs}}𝒢),$$
(8)
where $`𝒞`$ is the covariance matrix. Since $`\overline{\xi }^{\mathrm{mod}}_m`$ is a linear combination of template functions, the linear least squares estimators are

$$\widehat{𝒢}=\left(^T𝒞^1\right)^1^T𝒞^1\overline{\xi }^{\mathrm{obs}}$$
(9)
After obtaining $`𝒞`$, combined with an observed correlation function $`\overline{\xi }^{\mathrm{obs}}`$ and a design matrix $``$, we have the necessary machinery to obtain $`\widehat{𝒢}`$ i.e. the projection of $`\overline{\xi }^{\mathrm{obs}}`$ into the template functions. This means that separate fits for the lensing and intrinsic contributions can be obtained.

### 2.2 Covariance matrix

To evaluate (9), we need the covariance matrix $`𝒞_{mm^{}}`$. The covariance matrix could be calculated using the method described in Schneider et al. (2002), where $`𝒞`$ is expressed as integrals over (products of) correlation functions. However, since the method presented here would require the calculation of very many elements of the covariance matrix, owing to the redshift slicing, we decided to use, as a first step, a simplified model for $`𝒞`$. This consists of neglecting the cosmic variance contribution to $`𝒞`$, and thus consideration of the (diagonal) elements of $`𝒞`$ coming from the intrinsic ellipticity dispersion of the source galaxies. This is the dominant contribution to the covariance at small angular scales; at larger angular scales, the cosmic variance terms start to dominate, with the transition angular scale depending on the survey geometry (Kaiser 1998; Schneider et al. 2002). Here, we consider $`𝒩_{}`$ independent fields, and take $`𝒩_{}=300`$. Therefore, we expect the cosmic variance not to be very much larger than the intrinsic ellipticity noise on the angular scales considered.

With this approximation, the elements of the covariance matrix are

$$𝒞_{mm^{}}=\frac{2\left(\sigma _{ϵ}^{}{}_{}{}^{2}/2\right)^2}{𝒩_\mathrm{p}(m)}\left[\delta _{m,m^{}}\left(1+\delta _{I,I^{}}\delta _{J,J^{}}\delta _{I,J}\right)\right],$$
(10)
where there is an extra contribution from auto-variance terms. The elements of the inverse covariance matrix in this case are simply $`\left(𝒞^1\right)_{mm^{}}=\delta _{mm^{}}/𝒞_{mm}`$. The galaxy ellipticity dispersion is denoted by $`\sigma _ϵ`$. In bin $`m`$, the number of pairs is given by

$$𝒩_\mathrm{p}(m)=𝒩_\mathrm{F}\left[n_0p(z_I)\mathrm{\Delta }z\right]\left[n_0p(z_J)\mathrm{\Delta }z\right]L^4\frac{\mathrm{\Delta }\theta }{L}\tau \left(\frac{\theta }{L}\right),$$
(11)
where $`n_0`$ is the galaxy number density, $`p(z_I)`$ is the redshift probability density for redshift bin $`I`$, and $`L`$ is the extent of the field, assumed to be square. $`\tau (\theta /L)`$ is a function that takes into account the fact that fields have finite extent, and it must be evaluated numerically. Fig. 1 shows the function $`\tau (\theta /L)`$; note the limiting case where pairs separated by more than $`\sqrt{2}L`$ do not occur.

### 2.3 Basis models

In order to illustrate the general method described above, we will choose a simple, restricted set of template functions which share the approximate functional behaviour expected from the real correlation functions, both intrinsic and lensing. For the latter, we simply take a small number of CDM cosmologies and consider their correlation functions as a template set. For the former, simple exponentials with redshift dependence are chosen.

The lensing template functions $`A_n`$ used here are the gridded $`\overline{\xi }^\mathrm{L}(\theta ,\overline{z}_i,\overline{z}_j)`$ for 3 models of the underlying cosmology: $`\mathrm{\Lambda }`$CDM ($`\mathrm{\Omega }_\mathrm{m}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$), OCDM ($`\mathrm{\Omega }_\mathrm{m}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$) and $`\sigma `$CDM ($`\mathrm{\Omega }_\mathrm{m}=1`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$); for simplicity, $`\sigma _8=0.9`$ and $`\mathrm{\Gamma }=0.21`$ for each model. We use the Bardeen et al. (1986) transfer function to describe the evolution of the 3-dimensional power spectrum, along with the prescription of Peacock & Dodds (1996) for evolution in the non-linear regime. The required lensing correlation functions are calculated using the relationship between the power spectrum and $`\xi ^\mathrm{L}`$ given in (2), and then integrated over the photometric redshift uncertainties as in (4). Here it is assumed that $`p(z|\overline{z})`$ is a Gaussian with dispersion $`\sigma _{\mathrm{phot}}`$, centred on $`\overline{z}`$.

Nine template models $`B_n`$ for the intrinsic alignments are considered. First, the true spatial intrinsic correlation function is parameterised in terms of a correlation length $`R_{\mathrm{corr}}`$ and an exponent $`\alpha `$:

$$\eta (r,z)=(1+z_{\mathrm{av}})^\alpha \left[\mathrm{exp}\left(r/R_{\mathrm{corr}}\right)\right],$$
(12)
where $`z_{\mathrm{av}}`$ is the mean redshift of galaxies in a pair and $`r`$ is their comoving separation. We use the approximation $`r^2=(w_iw_j)^2+\theta ^2f^2\left[(w_i+w_j)/2\right]`$. $`R_{\mathrm{corr}}`$ was taken to be $`h^1\mathrm{Mpc}`$ and $`\alpha `$ to be \[$``$1, 0, 1\]. The availability of photometric redshift estimates is then accounted for by integrating $`\eta (r,z)`$ as in (4), and finally we obtain each of the model correlation functions on a grid. Note that the intrinsic models are calculated using the relationships for the distances $`f(w)`$ pertaining to the $`\mathrm{\Lambda }`$CDM cosmology described above, and that we use $`\xi _{\mathrm{gg}}(r)=(r/5h^1\mathrm{Mpc})^{1.8}`$.

## 3 Results

The results are presented in the context of a possible future multi-colour cosmic shear survey, with a field size $`L=14\mathrm{}`$, and $`𝒩_{}=300`$ independent pointings (i.e. the largest scale on which the ellipticity correlation function is available is $`\sqrt{2}\times 14\mathrm{}`$). A galaxy number density of 30 arcmin-2 and ellipticity dispersion of $`\sigma _ϵ=0.3`$ are used throughout. The value of $`\sigma _{\mathrm{phot}}=0.1`$ is chosen, typical of that obtained with current SED fitting procedures such as hyper-Z using a wide range of optical and near-infrared filters (Bolzonella et al. 2000). There are $`N_Z=65`$ redshift slices between $`\overline{z}=0.2`$ and 2.12, and $`N_K=25`$ angular separation bins between $`0.^{}3`$ and $`15\mathrm{}`$.

The galaxy redshift distribution follows the parameterisation suggested by Smail et al. (1995), i.e. $`p(\overline{z})=\beta /[z_0\mathrm{\Gamma }_\beta (3/\beta )](\overline{z}/z_0)^2\mathrm{exp}\left((\overline{z}/z_0)^\beta \right)`$, where $`\mathrm{\Gamma }_\beta `$ denotes the gamma function. We take $`\beta =3/2`$ and $`z_0=2/3`$ yielding $`\overline{z}1`$.

Two applications of the technique are presented. We start by asking how well the intrinsic and lensing contributions to the galaxy ellipticity correlation function can be separated, using the information contained in correlation functions between different redshift slices. Secondly, it has been shown that the degeneracy between $`\mathrm{\Omega }_\mathrm{m}`$ and $`\sigma _8`$ can be partly lifted when redshift estimates for source galaxies are available (e.g. van Waerbeke et al. 2002a). We illustrate the use of correlation function tomography in this respect.

### 3.1 Isolating the intrinsic correlation signal

Now we consider input “observed” correlation functions $`\overline{\xi }_m^{\mathrm{obs}}`$ comprising a lensing and an intrinsic contribution, to see how well the individual signals are recovered in terms of template functions using (9). First it was checked that in the absence of noise, $`\widehat{𝒢}_n𝒢_n`$, when $`\overline{\xi }^{\mathrm{obs}}`$ is composed of a lensing and an intrinsic model contained in the set of templates.

The intrinsic alignment model for spirals from HRH, $`\eta (r,z)=0.012\mathrm{exp}(r/1.5h^1\mathrm{Mpc})`$ was then used to obtain $`\overline{\xi ^\mathrm{I}}`$, and $`\overline{\xi ^\mathrm{L}}`$ was calculated for a $`\mathrm{\Lambda }`$CDM cosmology. Random gaussian distributed errors with dispersion $`\sigma =𝒞_{mm}^{0.5}`$ were added to these correlation functions giving noise realisations, and best-fit parameters $`𝒢_n`$ were recovered for each of these. Fig. 2 shows the (noise-free) input and recovered intrinsic and lensing correlation functions between three combinations of redshift slices: one close pair at $`\overline{z}1`$, one close pair at $`\overline{z}0.6`$, and the correlations between slices at $`\overline{z}1`$ and $`\overline{z}0.6`$. The lensing correlation function for each of the three cosmologies used in the construction of template functions is shown in Fig. 3, plotted for the slices at $`\overline{z}1.0`$. The intrinsic correlation signal surpasses the lensing signal out to several arcminutes for both the low-redshift and the high-redshift bins. Considering bins with a large separation in redshift reduces the intrinsic signal to a negligible level, as expected. Even with our limited set of template functions, the reduced $`\chi ^2`$ values of the recovered fits to the noise realisations are $`1`$. Also note that the intrinsic signal can be well represented in terms of the template functions, although it is not contained in the template set.

Since the current template set for the intrinsic alignment signal contains exponentially-decaying models with different scale-lengths, we now consider how well the method fares if the true signal is a power-law instead. The intrinsic alignment model is taken from Jing (2002); we use $`\eta (r,z)=0.288/[r^{0.4}(7.5^{1.7}+r^{1.7})]`$, with $`r`$ measured in units of $`h^1\mathrm{Mpc}`$. Again, the true cosmology is $`\mathrm{\Lambda }`$CDM. Noise was added and best-fit parameters recovered in the same manner as described above. Fig. 4 shows the (noise-free) input and recovered intrinsic and lensing correlation functions between the same three combinations of redshift slices as for Fig. 2. Even though the functional form of the true intrinsic signal is quite different from the template models, the best-fit intrinsic models are still rather close to the noise-free model and more importantly, the lensing signal is again well recovered. To assess the difference between using the Jing model rather than the HRH model for intrinsic alignments, the reduced $`\chi ^2`$ values of the recovered fits were determined for 1000 noise realisations of each, keeping the same $`\mathrm{\Lambda }`$CDM cosmology. The reduced $`\chi ^2`$ value is lower for the HRH realisation in nearly all cases, since this is an exponentially-decaying model for which the templates are better adapted. The $`\chi ^2`$ values for the HRH realisations closely follow the theoretically expected distribution. Although the intrinsic models can be distinguished statistically, the difference in the mean $`\chi ^2`$ values of the two sets of 1000 realisations is only $`15\%`$ of their dispersion. In practice, several families of functional forms could be taken for the template set.

The next example for $`\overline{\xi }^{\mathrm{obs}}`$ again uses the HRH model for intrinsic alignments, but this time a different flat cosmology with $`\sigma _8=0.71`$, $`\mathrm{\Omega }_\mathrm{m}=0.33`$, $`\mathrm{\Gamma }=0.215`$ was used for the lensing correlations; the lensing signal for this cosmology was not part of the template set. Fig. 3 shows the lensing correlation function for this cosmology, plotted for the the slices at $`\overline{z}1.0`$. Noise was added using the same random seeds as above, and best-fit parameters recovered as before. Fig. 5 shows results for the same combination of redshift slices as for Fig. 2. The power spectrum corresponding to this cosmology has a lower normalisation – the $`20\%`$ reduction in $`\sigma _8`$ is not offset by the $`10\%`$ increase in $`\mathrm{\Omega }_\mathrm{m}`$; this is evident in the lower amplitude of the lensing correlation functions. Again, the recovered parameters are consistent with the true “observed” correlation functions and the reduced $`\chi ^2`$ values $`1`$.

Whereas the foregoing cosmological model was quite different in amplitude of the power spectrum we now consider a model which also differs in the shape of the power spectrum, to test the robustness of our method. Keeping the HRH model for intrinsic alignments, a flat cosmology with $`\sigma _8=0.8`$, $`\mathrm{\Omega }_\mathrm{m}=0.5`$, $`\mathrm{\Gamma }=0.3`$ was next used for the lensing correlations. Fig. 3 shows the lensing correlation function for this cosmology, plotted for the the slices at $`\overline{z}1.0`$. Noise was added and best-fit parameters recovered as described above. Fig. 6 shows results for the same combination of redshift slices as for Fig. 2. The lensing signal is again well represented in terms of the basis functions, even though it is rather different to any of the cosmology templates. Hence, despite the fact that our set of template functions is quite restrictive, we have demonstrated that it provides enough flexibility to provide accurate fits to the correlation functions of quite different comological models.

### 3.2 Breaking the $`\mathrm{\Omega }_\mathrm{m}`$-$`\sigma _8`$ degeneracy

The near degeneracy, in the absence of redshift information, between $`\mathrm{\Omega }_\mathrm{m}`$ and $`\sigma _8`$ for two example flat cosmologies is illustrated in Fig. 7, taken from King & Schneider (2002; where details for its calculation can be found). One cosmology is the fiducial $`\mathrm{\Lambda }`$CDM model: $`\mathrm{\Omega }_\mathrm{m}=0.3`$, $`\sigma _8=0.9`$, and the other is an almost degenerate model with $`\mathrm{\Omega }_\mathrm{m}=0.4`$ and $`\sigma _8=0.78`$. To obtain these correlation functions, we assume the same prescription for the power spectrum outlined in Sect. 2.3 above. The source population has a redshift probability distribution with $`z=1`$, but with no individual photometric redshift estimates assumed to be available.

In order to see how well the fiducial and degenerate models could be distinguished using correlation function tomography, we took the gridded correlation functions for each model and added each of these to the gridded HRH intrinsic correlation function, giving $`\overline{\xi }^{\mathrm{mod}}_m`$ for each. For simplicity, we refer to these combinations as $`\overline{\xi }_{\mathrm{\Lambda }\mathrm{CDM}}`$ and $`\overline{\xi }_{\mathrm{degen}}`$. A first set of simulations involved using a set of ten template functions containing the nine models for intrinsic alignments, along with the $`\mathrm{\Lambda }`$CDM lensing model for the lensing template. In turn, 1000 noise realisations of $`\overline{\xi }_{\mathrm{\Lambda }\mathrm{CDM}}`$ and $`\overline{\xi }_{\mathrm{degen}}`$ were generated using the same random seeds in both cases, and the best-fit amplitudes $`\widehat{𝒢}_n`$ for the template functions recovered. The process was repeated, this time using the gridded degenerate model as the lensing template function, in place of the fiducial $`\mathrm{\Lambda }`$CDM model. The histograms of (i) $`\mathrm{\Delta }\chi ^2(\mathrm{\Lambda }\mathrm{CDM}\mathrm{degen})`$ and (ii) $`\mathrm{\Delta }\chi ^2(\mathrm{degen}\mathrm{\Lambda }\mathrm{CDM})`$, corresponding to the difference in goodness-of-fit for the noise realisations when $`\overline{\xi }_{\mathrm{\Lambda }\mathrm{CDM}}`$ and $`\overline{\xi }_{\mathrm{degen}}`$ are used in the template set, are shown in Fig. 8. When the fiducial $`\mathrm{\Lambda }`$CDM model (degenerate model) is the best-fit and is contained in the template set, values of $`\mathrm{\Delta }\chi ^2`$ in Fig. 8 should be negative. This gives a measure of our ability to differentiate between models using correlation functions between redshift slices. In the first histogram, when the model for $`\overline{\xi }_{\mathrm{\Lambda }\mathrm{CDM}}`$ is contained in the template set, in 95.6% of cases the noisy $`\mathrm{\Lambda }`$CDM correlation functions are better fit. Also, when $`\overline{\xi }_{\mathrm{degen}}`$ is in the template set, 96% of the noisy degenerate correlation functions have better fits. Hence, within the assumptions we made and in the presence of an intrinsic alignment signal, these two cosmological models could be distinguished at the $`2\sigma `$-level.

The foregoing example for distinguishing between two cosmological models should serve as a general illustration only. Owing to our simplified ansatz for the covariance matrix C, which contains intrinsic ellipticity dispersion only, a more detailed investigation is not warranted here. In practice, cosmic variance would need to be taken into account in the covariance matrix when realistic constraints on cosmological parameters are to be derived. We are currently investigating ways to obtain a far more realistic representation of the covariance matrix, to be used in a study of the accuracy of cosmological parameter determination.

## 4 Discussion and conclusions

It has been suggested that the lensing correlation function may be contaminated by intrinsic galaxy alignments. Since cosmic shear probes the matter power spectrum and enables constraints to be placed on cosmological parameters such as $`\sigma _8`$ and $`\mathrm{\Omega }_\mathrm{m}`$ (e.g. van Waerbeke et al. 2002a), it is vital to have the ability to isolate the contribution from intrinsic galaxy alignments in order to remove this systematic. Of course, intrinsic alignment is interesting in its own right: its amplitude as a function of physical separation and its evolution with redshift provides clues about the galaxy formation process.

We have demonstrated that measuring galaxy ellipticity correlation functions between redshift slices would enable the intrinsic and lensing contributions to be disentangled. The total signal is decomposed into template functions, and the fact that intrinsic alignments operate over a limited physical separation enables the intrinsic component to be isolated and subracted from the total signal. Our knowledge of the amplitude of intrinsic alignments is limited, but no strong assumption about the behaviour of the intrinsic alignment signal needs to be made. Here we considered a modest number of template functions, which can easily be augmented to cover a wider range of functional forms. For example, any intrinsic alignment signal arising at the epoch of galaxy formation may be suppressed by subsequent dynamical interaction, perhaps most pertinent to galaxy pairs with extremely small physical separations. In fact, if the reduced $`\chi ^2`$ of the best fit is significantly larger than 1, this indicates that additional template functions need to be included.

Our choice of template functions is of course fairly arbitrary. We have taken functions which have approximately the behaviour expected from a cosmic shear measurement. Alternatively, one could consider a set of generic basis functions, which however, owing to the dependence on three variables, would require a fairly large set of functions. Another natural choice of the template functions could be the following: assuming a reasonable guess for the cosmological model, characterised by the parameters $`\pi _0`$, the correlation function for neighbouring models could be written as

$`\overline{\xi }(\theta ,\overline{z}_i,\overline{z}_j;\pi )`$ (13)
$`\overline{\xi }(\theta ,\overline{z}_i,\overline{z}_j;\pi _0)+(\pi \pi _0){\displaystyle \frac{\overline{\xi }}{\pi }}(\theta ,\overline{z}_i,\overline{z}_j,\pi _0),`$
and therefore, the set consisting of $`\overline{\xi }(\pi _0)`$ and its partial derivatives with respect to the relevant cosmological parameters would provide a useful set of template functions. As an illustration, in Fig. 9 we have plotted the partial derivatives of $`\xi (\theta )`$ for the case where no redshift information is available, i.e. the derivatives of the redshift-averaged correlation function, again with $`\overline{z}1`$. For our cosmological model (characterised by $`\pi _0`$) we take the fiducial $`\mathrm{\Lambda }`$CDM model. Derivatives are taken with respect to (i) $`\mathrm{\Omega }_\mathrm{m}`$, (ii) $`\sigma _8`$, (iii) $`\mathrm{\Gamma }`$ and (iv) the primordial spectral index $`n`$ (where our fiducial model is scale-invariant $`n=1`$). We plot the ratio $`\frac{\xi }{\pi _i}/\xi `$, where the numerator is denoted by $`\mathrm{\Delta }(\xi ,\pi _i)`$. In the limiting case where one curve is a scaled version of another, $`\mathrm{\Delta }(\xi ,\pi _i)\mathrm{\Delta }(\xi ,\pi _j)`$, it is impossible to first order to break the degeneracy between parameters $`\pi _i`$ and $`\pi _j`$. We again see a nice illustration of the near degeneracy between $`\mathrm{\Omega }_\mathrm{m}`$ and $`\sigma _8`$, manifest in the similarity between $`\mathrm{\Delta }(\xi ,\mathrm{\Omega }_\mathrm{m})`$ and $`\mathrm{\Delta }(\xi ,\sigma _8)`$. Since $`\sigma _8^2`$ enters into the linear power spectrum just as a prefactor, on large angular scales (i.e. in the linear regime), the curve for $`\sigma _8`$ tends to a constant, $`\mathrm{\Delta }(\xi ,\sigma _8)/\xi 2\sigma _8^1`$ for large separations. Another feature to note is that on the scale of a few arcminutes the curves for $`\mathrm{\Gamma }`$ and $`n`$ change sign, implying that there is less degeneracy between these parameters and either of $`\mathrm{\Omega }_\mathrm{m}`$ or $`\sigma _8`$.

There are several ways in which the method and results discussed here could be used. One way would be to consider the resulting split into intrinsic correlations and lensing signal as the final result, and to compare the resulting functions with theories of galaxies formation which predict the intrinsic alignment signal, and cosmological models predicting the shear correlation function. The resulting fits are, however, difficult to interpret statistically, i.e. the error bars on the shear correlation function are difficult to obtain. An alternative would be to consider the fitted intrinsic signal only, subtract it from the ellipticity correlation function, and consider the result as the shear correlation function, together with the corresponding error bars. Subsequently, the correlation function can then be used for the redshift-weighting method of King & Schneider (2002), of course yielding much smaller contributions from the intrinsic correlations than for the unsubtracted data. Furthermore, the resulting model for the intrinsic correlation function could also be used as input for the subtraction method discussed in Heymans & Heavens (2002).

In addition to providing a key to the suppression of any intrinsic alignment signal, photometric redshift estimates enable much tighter constraints to be placed on cosmological parameters obtained from cosmic shear surveys, as demonstrated by Hu (1999). Although our prime goal in this paper is not the constraint of cosmological parameters, we have illustrated that the degeneracy between $`\mathrm{\Omega }_\mathrm{m}`$ and $`\sigma _8`$ can be lifted by observing correlation functions between redshift slices, even when an intrinsic alignment systematic is present.

###### Acknowledgements.

We would like to thank Marco Lombardi, Patrick Simon, Douglas Clowe, and in particular Martin Kilbinger for helpful discussions. We would also like to thank the anonymous referee for very helpful remarks. This work was supported by the Deutsche Forschungsgemeinschaft under the project SCHN 342/3–1, and by the German Ministry for Science and Education (BMBF) through the DLR under the project 50 OR 0106..