text:
# Spontaneous Vortex Phase in the Bosonic RVB Theory

## I Introduction

A unique feature of superconductivity in the high-$`T_c`$ cuprates is that phase coherence in electron pairing may occur at a lower temperature than the temperature at which the pairing amplitude develops . Doped Mott insulators provide a natural explanation: owing to the separation of spin and charge degrees of freedom, the spin resonating valence bond (RVB) pairing can be achieved at a much higher temperature; while the RVB pairs, accompanied by charge backflow, can move around like Cooper pairs at finite doping, superconducting condensation is absent until phase coherence between the pairs is established at a relatively lower $`T_c`$ .

Recently, we proposed a Ginzburg-Landau description of the RVB superconductor  based on the phase string theory of the $`tJ`$ model. In this description, the superconducting order parameter is given by

$$\mathrm{\Delta }(𝐫)\mathrm{\Delta }^0(𝐫)e^{i\mathrm{\Phi }^s(𝐫)}$$
(1)
where the Cooper-pair amplitude $`\mathrm{\Delta }^0\mathrm{\Delta }^s[\psi _h^{}]^2,`$ with $`\mathrm{\Delta }^s`$ being the bosonic RVB order parameter of spins and $`\psi _h`$ being the charge (holon) Bose condensed field. Here $`\mathrm{\Delta }^00`$ does not directly result in superconductivity; rather, it is the phase factor in (1) that determines the phase coherence of the pairing order parameter, and thereby, $`T_c`$.

The quantity $`\mathrm{\Phi }^s(𝐫)`$ characterizes phase vortices centered around spinons, since $`\mathrm{\Phi }^s(𝐫)\mathrm{\Phi }^s(𝐫)\pm 2\pi `$, if $`𝐫`$ winds around a spinon excitation continuously in space. Thus, each spinon excitation induces a phase vortex (called spinon vortex) in the order parameter. Superconducting phase coherence is destroyed above $`T_c`$ by the presence of free spinon vortices. Below $`T_c`$, phase coherence is realized as spinon vortices and antivortices are bound , such that $`\mathrm{\Phi }^s(𝐫)`$ becomes trivial in (1). In the superconducting phase, single spinon vortices can only be present at magnetic vortex cores which ensures flux quantization at $`hc/2e`$.

It should be noted that the bosonic RVB order parameter $`\mathrm{\Delta }^s`$ in (1) is not related to the energy gap, in contrast to the usual fermionic RVB order parameter. It describes (neutral) spin pairing as characterized by short range (nearest-neighbor) antiferromagnetic correlations, $`𝐒_i𝐒_j_{NN}=1/2|\mathrm{\Delta }^s|^2`$. At small doping, $`\mathrm{\Delta }^s0`$ covers a temperature regime extended over $`1,000`$K . On such an RVB background, the Cooper pair amplitude $`\mathrm{\Delta }_0=\mathrm{\Delta }^s[\psi _h^{}]^2`$ is realized when the holons (Bose) condense, i.e., $`\psi _h0`$. But the temperature $`T_v`$ at which this occurs does not necessarily coincide with $`T_c`$ in general. Thus, one may find a temperature regime in the normal state, $`T_v>T>T_c`$ (obviously $`T_v`$ cannot be lower than $`T_c`$ since it is the holon Bose condensation that engenders phase coherence ).

In this paper, we will focus on a normal state with nonzero (preformed) pairing amplitude. It is distinguished from a conventional normal state of strong superconducting fluctuations by the presence of free spinon vortices as elementary excitations. We call such a phase, a spontaneous vortex phase. We present the effective Hamiltonian governing the dynamics of spinon vortices. Some novel transport properties in this state are discussed, and in particular, we show that free spinon vortices contribute to a nontrivial Nernst signal, consistent with the recent measurements of the Princeton group . We argue that a finite $`\mathrm{\Delta }_0`$ controls the pseudogap phenomena, and present a phase diagram showing the connection between the superconducting phase and the phase with preformed pairs (spontaneous vortex phase).

## II Spinon Vortices: Elementary Excitations

### A Spinons as vortices

The phase $`\mathrm{\Phi }^s(𝐫)`$ in (1) is given by

$$\mathrm{\Phi }^s(𝐫)=d^2𝐫^{}\text{Im}\mathrm{ln}\left[zz^{}\right]\left[n_{}^b(𝐫^{})n_{}^b(𝐫^{})\right]$$
(2)
where $`n_\sigma ^b`$ denotes the density operator of bosonic spin-1/2 excitations (spinons) and $`z=x+iy`$. As noted in the Introduction, $`\mathrm{\Phi }^s(𝐫)`$ changes by $`\pm 2\pi `$ when $`𝐫`$ winds around a spinon once. In the ground state, the spinons are all RVB paired such that $`\mathrm{\Phi }^s(𝐫)`$ is effectively canceled out. Excited spinons created by breaking the RVB pairs lead to vortices in the pairing order parameter (1) through $`\mathrm{\Phi }^s(𝐫)`$.

The spinon vortex can be understood from a different perspective. The effective Hamiltonian for the holons is given by

$$H_h=\frac{1}{2m_h}d^2𝐫h^{}(𝐫)\left(i𝐀^s𝐀^e\right)^2h(𝐫)$$
(3)
where $`h(𝐫)`$ is the bosonic holon field, and $`𝐀^e`$ is the vector potential of the external electromagnetic field. (Here we set $`\mathrm{}=e=c=1`$. ) The quantity $`𝐀^s`$ defined by

$$𝐀^s(𝐫)=\frac{1}{2}d^2𝐫^{}\frac{\widehat{𝐳}\times (𝐫𝐫^{})}{|𝐫𝐫^{}|^2}\left[n_{}^b(𝐫^{})n_{}^b(𝐫^{})\right],$$
(4)
is related to $`\mathrm{\Phi }^s`$ by $`𝐀^s=\frac{1}{2}\mathrm{\Phi }^s`$. Note that

$$_c𝑑𝐫𝐀^s=\frac{1}{2}_c𝑑𝐫𝚽^s=\pm \pi _{Sc}d^2𝐫^{}\left[n_{}^b(𝐫^{})n_{}^b(𝐫^{})\right]$$
(5)
i.e., $`𝐀^s`$describes a fictitious flux quantized at $`\pm \pi `$, bound to each spinon as seen by holons.

The vortices in $`\mathrm{\Phi }^s(𝐫)`$ become meaningful when $`\mathrm{\Delta }^00`$ in (1), i.e., when the holons are Bose condensed. In terms of the superfluid density $`\rho _h`$ and the phase $`\varphi _h`$, the holon condensate $`\psi _h(𝐫)h(𝐫)=\sqrt{\rho _h}e^{i\varphi _h(𝐫)}`$, and equation (3) reduces to

$$H_h\frac{1}{2m_h}d^2𝐫\rho _h\left(\varphi _h𝐀^s𝐀^e\right)^2$$
(6)
The corresponding superfluid current operator is given by

$$𝐉=\frac{\rho _h}{m_h}\left[\varphi _h𝐀^s𝐀^e\right].$$
(7)
Consider a spinon excited from the RVB background within a loop $`c`$. In terms of (7), there must be an induced current vortex surrounding the excited spinon (setting $`\varphi _h=𝐀^e=0)`$

$$_c𝐉(𝐫)𝑑𝐫=\frac{\rho _h}{m_h}_c𝐀^s(𝐫)𝑑𝐫=\frac{\rho _h}{m_h}\pi .$$
(8)
Therefore, each spinon excitation will always be accompanied by a current vortex in the holon condensate. This current vortex is consistent with the aforementioned $`\pm 2\pi `$ vorticities of $`\mathrm{\Phi }^s`$ in $`\mathrm{\Delta }`$ (in the latter a pair of holon fields are present).

It should be noted that such a vortex cannot be compensated or screened by $`\varphi _h`$ in (7) since the latter must satisfy the single-valued condition $`𝑑𝐫\varphi _h=0`$, $`\pm 2\pi `$, …, instead of $`\pm \pi `$, the flux quanta carried by spinons in $`𝐀^s`$. However, $`\varphi _h`$ provides a gauge freedom such that the sign of the vorticity for each minimal current vortex surrounding an excited spinon in (8) is not related to the spin polarization itself . For example, consider a spinon at site $`i`$, with $`𝑑𝐫𝐀^s=\frac{1}{2}𝑑𝐫𝚽^s=\pi `$. Under a singular gauge transformation, $`\psi _h(𝐫)\psi _h(𝐫)e^{i\theta (𝐫)}`$ with $`𝑑𝐫\theta =2\pi `$ centered at $`i`$, the theory is invariant if $`𝐀^s𝐀^s\theta `$ ($`\mathrm{\Phi }^s\mathrm{\Phi }^s2\theta `$), with $`𝑑𝐫𝐀^s=\frac{1}{2}_c𝑑𝐫\mathrm{\Phi }^s\pi `$. In the appendix, we give a general proof that spin rotational symmetry is indeed preserved here.

Finally, we note that there also exist the usual Kosterlitz-Thouless (KT) like $`2\pi `$ vortices in $`\psi _h.`$ But since we will only be interested in the temperature regime, $`T<T_{v\text{ }}`$ where holons are condensed, we ignore their effects. Consequently, in the rest of the paper, $`\varphi _h`$ will be set to zero in (6) and (7), with its singular part associated with spinons being incorporated into $`𝐀^s`$. So, in the latter \[(4)\] the spin index should be generally construed as the vorticity index and independent of the actual spin index, as discussed above.

### B Effective Hamiltonian of spinon vortices

The effective Hamiltonian governing the spinon vortices can be written in two parts:

$$H_{spinonvortex}=H_s+H_v$$
(9)
Here, $`H_s`$ is the mean field Hamiltonian for bosonic spinon excitations,

$$H_s=\underset{m\sigma }{}E_m\gamma _{m\sigma }^{}\gamma _{m\sigma }+const.$$
(10)
where $`\gamma _{m\sigma }`$ is related to the bare bosonic operator by the Bogoliubov transformation

$$b_{i\sigma }=\underset{m}{}w_{m\sigma }\left(𝐫_i\right)\left[u_m\gamma _{m\sigma }v_m\gamma _{m\sigma }^{}\right]$$
(11)
with $`u_m=1/\sqrt{2}(\lambda _m/E_m+1)^{1/2}`$ and $`v_m=1/\sqrt{2}(\lambda _m/E_m1)^{1/2}sgn(\xi _m)`$; $`E_m=\sqrt{\lambda _m^2\xi _m^2}`$  and $`\lambda _m=\lambda \frac{J_h}{J_s}\left|\xi _m\right|`$ ($`J_h\delta t`$, $`J_sJ)`$. The single particle wave function $`w_{m\sigma }\left(𝐫_i\right)`$ and $`\xi _m`$ are determined by the tight binding equation

$$\xi _mw_{m\sigma }\left(𝐫_i\right)=J_s\underset{j=NN(i)}{}e^{i\sigma A_{ji}^h}w_{m\sigma }(𝐫_j)$$
(12)
in which $`w_{m\sigma }\left(𝐫_i\right)=w_{m\sigma }^{}\left(𝐫_i\right).`$ Note that $`_cA_{ij}^h=\pi _{lc}n_l^h`$ describes $`\pi `$ flux tubes bound to holons ($`n_l^h`$ denotes the holon number operator). In the holon Bose condensed phase, it is a good approximation to treat $`A_{ij}^h`$ as a vector potential for a uniform flux perpendicular to the two-dimensional plane, with a field strength

$$B^h=\frac{\pi \delta }{a^2}$$
(13)
where $`\delta `$ is the doping concentration and $`a`$ is the lattice constant.

Therefore, the tight binding model in (12) has a Hofstadter-Landau level structure. In the weak-field (continuum) limit, $`\xi _m`$ is a discrete, dispersionless Landau-level energy, with the wave function characterized by a cyclotron length scale

$$a_c=\frac{1}{\sqrt{B^h}}=\frac{a}{\sqrt{\pi \delta }}$$
(14)
Correspondingly the spinon spectrum $`E_m`$ is discrete as well, and if we restrict ourselves to the lowest Landau-level (LLL) at the low temperature, the dispersionless spectrum

$$(E_m)_{\text{LLL}}E_s$$
(15)
with a degeneracy equal to $`2\times B^ha^2/2\pi =`$ $`\delta `$ (the prefactor $`2`$ comes from the fact that $`\xi _m`$ and -$`\xi _m`$ solutions are degenerate in $`E_m)`$ and $`E_s\delta J`$  ($`\lambda `$ in $`E_m`$ is determined by the average constraint condition $`_\sigma b_{i\sigma }^{}b_{i\sigma }=1\delta `$). In this theory,

$$E_g=2E_s$$
(16)
corresponds to the energy scale of the sharp resonance-like peak observed in the neutron-scattering measurements, with its weight proportional to the degeneracy of the energy levels, $`\delta `$.

The second term in (9), $`H_v`$, describes the interaction among spinons due to the fact that they carry current vortices. This term arises from the effective Hamiltonian for the holons, $`H_h`$. Setting $`\varphi _h=𝐀^e=0`$ in (6), we get

$`H_v`$ $`=`$ $`{\displaystyle \frac{\rho _h}{2m_h}}{\displaystyle d^2𝐫\left(𝐀^s\right)^2}`$ (17)
$`=`$ $`{\displaystyle d^2𝐫_1d^2𝐫_2_\alpha \alpha n_\alpha ^b(𝐫_1)V(𝐫_{12})_\beta \beta n_\beta ^b(𝐫_2)}`$ (18)
in which

$`V(𝐫_{12})`$ $`=`$ $`{\displaystyle \frac{\rho _h}{2m_h}}{\displaystyle d^2𝐫\frac{\widehat{𝐳}\times (𝐫𝐫_1)}{|𝐫𝐫_1|^2}\frac{\widehat{𝐳}\times (𝐫𝐫_2)}{|𝐫𝐫_2|^2}}`$ (19)
$`=`$ $`{\displaystyle \frac{\pi \rho _h}{4m_h}}\mathrm{ln}{\displaystyle \frac{|𝐫_1𝐫_2|}{r_c}}`$ (20)
with $`r_ca`$. Using the Bogoliubov transformation (11) and defining $`n_{m\sigma }^\gamma =\gamma _{m\sigma }^{}\gamma _{m\sigma }`$ , one finds

$$_\sigma \sigma n_\sigma ^b(𝐫)=\underset{m}{}\left|w_{m\sigma }(𝐫)\right|^2_\sigma \sigma n_{m\sigma }^\gamma +\underset{mn}{}\mathrm{}$$
(21)
and thus the interaction term can be further rewritten as $`H_v=H_v^0+H_v^1,`$ with

$$H_v^0=\underset{\alpha \beta }{}\alpha \beta \underset{mn}{}n_{m\alpha }^\gamma U_{mn}n_{n\beta }^\gamma $$
(22)
where

$$U_{mn}=d^2𝐫_1d^2𝐫_2\left|w_{m\alpha }(𝐫_1)\right|^2\left|w_{n\beta }(𝐫_2)\right|^2V(𝐫_{12})$$
(23)
and the non-diagonal part $`H_v^1`$ describes the scattering among different states induced by vortex interactions.

The diagonal part $`H_v^0`$ provides a confining force for spinons in the superconducting phase. Unlike the conventional KT vortices, there is a finite core for each spinon-vortex composite in $`H_v^0`$, within which the spinon does a cyclotron motion with a core radius $`a_c`$ determined by $`\left|w_{m\alpha }(𝐫)\right|^2`$. From a detailed renormalization group analysis, it is found that $`T_c`$, the temperature at which spinon vortices unbind, scales with the spin resonance-like energy $`E_g`$, consistent with a simpler physical argument based on the “core touching” picture, where it has been estimated that the transition occurs when the number of excited spinons exceeds the holon number.

## III Spontaneous vortex phase

### A Definition

The spinon-vortex composite described above is a unique elementary excitation, predicted by the present bosonic RVB theory, in the regime $`\psi _h`$ or $`\mathrm{\Delta }_00`$ below the holon condensation temperature $`T_v.`$ As pointed out in the Introduction, this regime generally includes two phases. One is the superconducting phase at a lower temperature $`T_c`$, where phase coherence is realized when spinon vortices are paired up (“confined”). The second phase at $`T_c<T<T_v`$ is a normal state in which unbinding (“deconfined”) free spinon vortices are present. We define such a state of matter as the spontaneous vortex phase. One may also regard this phase as a normal state with a finite Cooper-pair amplitude $`\mathrm{\Delta }_0`$.

### B Phase rigidity and the Nernst effect

The existence of free vortices implies that some kind of phase rigidity persists above $`T_c`$ in the spontaneous vortex phase. In the following, we explore some consequences.

Consider a strip sample showing in Fig. 1, in which a spinon vortex traverses from one end of the sample to the other along the strip. According to (1) and (2), it is straightforward to see that the phase difference in $`\mathrm{\Delta }`$ between two edges across the strip will change by $`2\pi `$ (phase slip) for an infinitely long strip. Hence the boundaries of the sample can always perceive the motion of spinon vortices inside, which is a direct indication of the phase rigidity. The presence of such a new elementary excitation is the most important feature distinguishing the spontaneous vortex phase from an ordinary normal state with strong superconducting fluctuations.

A phase slip of $`2\pi `$ in a Josephson junction represents an elementary resistivity process as a voltage will be generated in that (transverse) direction. In an equilibrium state at $`T_c<T<T_v`$, phase slippage is random as those free spinon vortices move in arbitrary directions (this is an another way to understand the destruction of the phase coherence in this regime) and the average voltage is zero. In order to see a net phase slippage being added up, one has to let all the spinon vortices move in the same direction, say, by maintaining a temperature gradient and use a perpendicular magnetic field to induce an imbalance between vortices and antivortices. This results in the well known Nernst effect. The existence of a nontrivial contribution from spinon vortices to the Nernst signal, is one of the main features of the spontaneous vortex phase.

There is another way to perceive this effect. In the superconducting phase, it is well known that the motion of magnetic vortices generates the Nernst signal. In the present case, as discussed in Ref., a magnetic vortex is formed by a spinon bound to a magnetic flux quantized at $`hc/2e`$ (equal to $`\pi `$ in the present units$`)`$ \[Fig. 2(a)\]. This can be easily seen based on (1) which shows that each $`2\pi `$ vortex in the superconducting order parameter must be associated with a spinon through $`\mathrm{\Phi }^s.`$ In the bulk, the spinon vortices have to be paired and do not contribute to $`\mathrm{\Phi }^s`$ and the net Nernst effect. So the Nernst signal from $`\mathrm{\Phi }^s`$ arises only from those spinon vortices that are nucleated by the magnetic fluxoids and are bound to the latter as shown in Fig. 2(a).

At $`T_c,`$ the phase coherence is destroyed as $`\mathrm{\Phi }^s`$ becomes disordered in (1) due to the unbinding of spinon vortices, and $`e^{i\mathrm{\Phi }^s(𝐫)}e^{i\mathrm{\Phi }^s(𝐫^{})}`$ falls off exponentially at large distance. Correspondingly, the magnetic vortices collapse and the Meissner effect disappears; the spinons originally trapped to the magnetic vortex cores are released from the latter and become free, as illustrated in Fig. 2(b). Since, in the present theory, the spinon vortices contribute to the Nernst signal through the order parameter phase $`\mathrm{\Phi }^s`$, the collapse of magnetic vortices themselves does not directly lead to a diminishing Nernst signal above $`T_c`$. On the other hand, unbinding vortices and vortices do not give rise to additional  Nernst contribution since the net imbalance of vortices and antivortices

remains the same below and above $`T_c`$, determined by the condition

$$𝐉(𝐫)𝑑𝐫=\mathrm{𝟎}$$
(24)
or

$$𝑑𝐫𝐀^e=𝑑𝐫𝐀^s.$$
(25)
The only difference is that below $`T_c`$, the above condition holds for a closed loop encircling a magnetic fluxoid and far away from the core. (It also leads to flux quantization at $`hc/2e`$ as discussed in Ref..) Above $`T_c`$, without magnetic vortices, it holds on a spatial average. Thus the Nernst signal arises from the same spinons trapped at magnetic vortex cores below $`T_c`$ \[Fig. 2(a)\] and released in the spontaneous vortex phase \[Fig. 2(b)\], above $`T_c`$. Therefore, we expect the Nernst effect evolve smoothly between the superconducting and spontaneous vortex phases at $`TT_c.`$

This smooth evolution of the Nernst effect is probably the most striking feature observed in the recent Nernst measurement on LSCO compounds by Xu et al.. These measurements clearly reveal a continuous crossover of the Nernst signal above $`T_c`$; it remains anomalously enhanced up to temperatures, $`50100`$ K above $`T_c`$. Of course, experimentally one has to separate the vortex contribution to the Nersnt signal from that of normal charge carriers whose contribution is limited by the so called Sondheimer cancellation. Xu et al.  argue that the Nersnt signal arises from vortex like excitations in the normal state and are not likely to be related to the conventional superconducting fluctuations. Similar phenomena have also been found in YBCO and BSCO systems as well, indicating that the existence of vortices above $`T_c`$ is a generic feature of the high-$`T_c`$ cuprates.

### C Transport coefficients

#### 1 Spinon vortices

Spinons do not carry charge and thus do not couple to the external electromagnetic field $`𝐀^e`$ directly. However, they are vortices and thus generate a Nernst voltage in the direction transverse to their motion, due to the phase slip effect discussed above. This can be quantified as follows. In terms of (7), one has

$$_t𝐉=\frac{\rho _h}{m_h}\left(_t𝐀^s+_t𝐀^e\right).$$
(26)
The electric field is given by $`𝐄=_t𝐀^e`$ in the transverse gauge (with $`\varphi _h`$ being absorbed by $`𝐀^e,`$ which is always possible below $`T_v`$). In the steady state, $`_t𝐉=0`$, and we get

$$𝐄=\underset{l}{}\widehat{𝐳}\times 𝐯^s(l)\alpha _l\delta (𝐫𝐫_l)$$
(27)
where $`l`$ labels the spinon vortices with $`\alpha _l=\pm \pi `$ denoting the vorticity. Equation (27) means that a finite drift velocity $`v_x^s`$ of a spinon-vortex along $`\widehat{x}`$-direction indeed will induce an electric voltage along the $`\widehat{y}`$-direction, as the result of the phase slip effect shown in Fig. 1. Using the condition (25), we get

$$𝐁=\widehat{𝐳}\underset{l}{}\alpha _l\delta (𝐫𝐫_l)$$
(28)
where $`𝐁`$ denotes the local magnetic field perpendicular to the 2D plane. If all vortices have the same drift velocity, we get $`𝐄=𝐁\times 𝐯^s`$ which coincides with the familiar form for magnetic fluxoids in the flux flow region of a type II superconductor.

The Nernst coefficient is defined as

$$\nu _{sv}=\frac{E_y}{(_xT)B_z}.$$
(29)
Note that in this measurement $`J_y`$ is set to be zero. $`\nu _{sv}`$ can be determined phenomenologically. Suppose $`s_\varphi `$ is the transport entropy carried by a spinon vortex and $`\eta _s`$ is its viscosity. Then the drift velocity $`𝐯^s`$ can be decided by

$$s_\varphi T=\eta _s𝐯^s,$$
(30)
and consequently,

$$\nu _{sv}=\frac{s_\varphi }{\eta _s}.$$
(31)
The main result of the bosonic RVB theory is that such an expression is meaningful both above and below $`T_c`$. As discussed earlier, $`\nu _{sv}`$ evolves smoothly between the superconducting and spontaneous vortex phase at $`TT_c`$, since the Nernst signal arises from the spinon vortices in (27), and the formation or collapse of magnetic vortices does not change the expression (31). Of course, below $`T_c`$, $`s_\varphi `$ will include additional contributions associated with the magnetic vortex core, which remains to be investigated.

Thus, the existence of spinon vortices provides a natural explanation for the Nernst experiments in the spontaneous vortex phase. In this case, the Nernst effect directly measures the dynamics of spinon vortices as governed by (9). The viscosity, $`\eta _s`$, of spinon vortices is proportional to the scattering rate of the spinons. The latter is enhanced above $`T_c`$ due to the broadening of the spin resonance level at $`E_g`$, and correlated with the fate of the resonance peak observed in inelastic neutron scattering measurements. As will be discussed in the next subsection, $`\eta _s`$ is also related to the resistivity.

The transport entropy $`s_\varphi `$, is related to the degeneracy of the level $`E_s`$ as the “mid-gap” states of the spinon trapped inside the magnetic vortex core. We expect it to decrease rapidly above $`T_c`$. This is not due to the collapse of the bound core states (of the trapped spinon), as in the conventional superconductors. It is due to the 2D Coulomb interaction \[(18)\], whose effect is important in the high (vortex) density limit. Note that the degeneracy of the spinon states at the level $`E_s`$ is related to different cyclotron orbits in the LLL (see Sec. IIB). In the “core touching” picture, the number of excited spinons at $`T_c`$, is approximately equal to the holon number such that each LLL state, corresponding to the resonance peak at $`E_g,`$ is occupied approximately by one spinon on average. Higher temperature only means more spinon vortices will be excited to the LLLs such that more than one spinon will be found within a cyclotron orbit. In this dense limit, those configurations with many vortices (of the same vorticity) lumped together, which are allowed by the bosonic statistics, are forbidden due to costing too much potential energy. Thus, the phase space of spinons will be strongly limited in the dense limit, as compared to the dilute case, due to the vortex-vortex interaction, and the entropy associated with each spinon vortex should be largely reduced above $`T_c`$ such that the Nernst coefficient $`\nu _{sv}`$ would decrease in magnitude rapidly. A microscopic study based on the effective Hamiltonian (9) is needed in order to understand this issue and make quantitaive comparison with experimental results.

Spinon vortices also contribute to the thermal current established by a temperature gradient,

$$𝐉^Q=\kappa _{sv}(T)$$
(32)
Here the heat current is given by

$$𝐉^Q=\underset{m\sigma }{}𝐯^s(m)E_mn_{m\sigma }^\gamma 𝐯^sE_sn_v$$
(33)
with $`n_v=_{m\sigma }n_{m\sigma }^\gamma `$. In obtaining the last step in (33), it has been assumed that spinon vortices are mainly excited to the lowest discrete energy level at $`E_s`$ and the interaction among vortices in (18) only gives rise to level broadening which is averaged out in (33) with $`𝐯^s`$ representing an average drift velocity.

Based on (30)-(33), the following connection between the thermal conductivity and Nernst effect of the spinon vortices can be deduced:

$$\frac{\kappa _{sv}}{\nu _{sv}}=\left(\frac{E_g}{2}\right)n_v.$$
(34)
At $`T_c`$, one has $`n_v\delta a^2`$ , and for optimally doped YBCO with $`E_g41`$meV and $`\delta 0.15`$, we estimate $`\kappa _{sv}/\nu _{sv}3`$meV$`/a^2`$.

Finally, as a prediction, we introduce

$$\zeta _{sv}J_x^Q/E_yB_z$$
(35)
in the following experiment setting: apply an electric field along $`\widehat{y}`$-direction and measure the heat current along $`\widehat{x}`$-direction in the presence of a perpendicular magnetic field $`B_z`$. In such an experimental situation, we expect that spinon vortices contribute predominantly to $`J_x^Q`$. No contribution from phonons will be generated by $`E_y`$, while the Hall effect and the thermal conductivity of the charge carriers, holons, are small in the spontaneous vortex phase as will be discussed shortly. In this case, an electric current induced by $`E_y`$ along the $`\widehat{y}`$-direction exerts a Lorentz force on the vortices causing them to move in the $`\pm \widehat{x}`$-directions, depending on their vorticities (see discussion in the next subsection). Then, (27) and (28) can be rewritten as $`E_y=v_x^s\pi n_v`$ and $`B_z=\pi \mathrm{\Delta }n_v`$, respectively, where $`\mathrm{\Delta }n_v`$ denotes the difference between the numbers of positive and negative vortices. The thermal current $`J_x^Q=v_x^sE_s\mathrm{\Delta }n_v`$, is related to $`\mathrm{\Delta }n_v`$ instead of $`n_v`$, because vortices with opposite vorticities move in opposite directions. Then, we obtain

$$\zeta _{sv}=\frac{2E_g}{\mathrm{\Phi }_0^2n_v},$$
(36)
where $`\mathrm{\Phi }_0/2hc/2e`$ (full units have been restored). This quantity is independent of the transport entropy and viscosity, and is only related to the characteristic spin excitation energy scale and the number of spinon vortices.

#### 2 Holons: Charge carriers

The bosonic holon carries charge $`+e`$ and directly couples to the electromagnetic field in (3). In the spontaneous vortex phase, holons are Bose-condensed. But the “superfluid” density does not carry zero resistivity in the spontaneous vortex phase. Rather, one expects to see a finite resistivity from the condensed holons, which arises from the motion of spinon vortices. Suppose an electric current is established along the $`\widehat{x}`$-direction. Each spinon vortex carries a $`\mathrm{\Phi }_0/2`$ fictitious flux with a radius $`a_c`$ in the holon Hamiltonian (3). Like a conventional magnetic vortex in a type II superconductor, a transverse “Lorentz force” due to the interaction of the current with the fluxoid deflects the spinon vortex along the $`\pm \widehat{y}`$ directions, depending on the sign of vorticity. In turn, the moving spinon vortex produces a voltage along the $`\widehat{x}`$ direction, as can be seen from (27). It is then easy to obtain

$$\rho =\frac{n_v}{\eta _s}\left(\frac{\mathrm{\Phi }_0}{2c}\right)^2.$$
(37)
On the other hand, due to the condition (25), the magnetic field seen by holons will be canceled out in (3) on the average by $`𝐀^s`$. Thus the Hall resistivity will be reduced in the spontaneous vortex phase. By the same token, the magnetoresistance should also be weak, as holons do not see a net field on an average. Note that $`\rho `$ is not sensitive to how vortices and antivortices get polarized by a magnetic field. As seen from (37), the resistvity depends only on the total number of excited spinon vortices, $`n_v,`$ which as a function of $`T`$ and $`E_s`$ should not be very sensitive to the magnetic field, at least in the weak field limit. Finally, due to the fact that the holon condensate does not carry an entropy, the thermopower is also expected to be suppressed below $`T_v`$ and only the normal-fluid component has a residual contribution.

Transport measurements of many cuprate superconductors show that the suppression of the Hall effect and thermopower around the same temperature scale ($`T_v`$) above $`T_c`$ where the vortex-induced Nernst effect ends. These observations are complemented by a weak magnetoresistance in the same regime. All of this lends a support for the existence of the spontaneous vortex phase below $`T_v`$, which can provide a consistent picture as discussed above.

So far, we have not discussed the role of the nodal (fermionic) quasiparticle. In the bosonic RVB theory, the quasiparticle arises as a composite object of a confined holon-spinon pair, with a $`d`$-wave dispersion  like in a $`d`$-wave BCS superconductor. But above $`T_c`$, due to the deconfinement of vortices and antivortices, such a quasiparticle is expected to be damped, as its spinon constituent can get away freely. Thus, except for residual effects, their contribution to transport should become negligible for $`T>T_c`$. In contrast, the contribution of the quasiparticles becomes dominant below $`T_c`$, when their coherence is established.

## IV Phase diagram and discussions

The contribution of spinon vortices to the Nernst effect vanishes at $`T_v`$, when holon condensation or $`\mathrm{\Delta }_0`$ disappears. In this section, we discuss the evolution of the spontaneous vortex phase as a function of doping and temperature and propose a phase diagram.

Physically, holon condensation at $`T_v`$ can be interpreted as the original neutral RVB pairs at high temperature acquiring charge and becoming Cooper-like pairs. Above $`T_v`$, incoherent holons are strongly scattered by the spinon flux-tubes according to (3). So even though the motion of an RVB spin pair may be accompanied a pair of holons as a backflow, it does not mean that the RVB pair is “charged” by -2e. The 2e holon backflow becomes coherent only below $`T_v`$, when the two holons constituting a pair are condensed. In this regime, the RVB spin pair looks like carrying a charge -2e (with a single wavefunction). This is the physical picture of how the amplitude of the pairing parameter forms below $`T_v`$.

The dashed line in Fig. 3 shows the bare KT temperature $`T_{KT}`$ for the holon condensation, which is obtained without including $`𝐀^s`$ in (3). First, we argue that $`T_vT_{KT}`$ at very low doping concentrations. To understand the effect of fluctuating $`𝐀^s,`$ let us first recall the physical interpretation of the $`T_c`$ curve in Fig. 3. The superconducting transition occurs due to the dissolution of spinon vortex-antivortex pairs in (1). Such a transition is not driven by the entropy reason as in a conventional KT transition. Rather, it is related to the so-called “core touching” mechanism: each spinon vortex has a finite core with the length scale compatible to $`a_c`$, and an upper limit of $`T_c`$ is determined by the temperature at which the cores of excited spinon-vortices start to touch, which leads to $`T_cE_g/4`$. Considering that the overlapping of vortex cores effectively smear out the flux fluctuations in (4), in spite of unbinding spinon vortices, the fluctuations of $`𝐀^s`$ are not necessarily strong above $`T_c`$. This is particularly true at small doping when $`a_c`$ can become very large. Hence, one expects that $`T_v`$ approximately coincides with the bare KT temperature, $`T_{KT}=\pi \delta (2a^2m_h)^1`$, for small doping.

With the increase of doping, the core scale gets reduced. Near optimal doping, $`\delta =0.15,`$ $`a_c`$ $`1.5a`$ which becomes comparable to the lattice constant. At such a short length scale, core overlapping no longer results in the smoothness of $`𝐀^s`$. Instead, the uncorrelated spinons above $`T_c`$ gives rise to strong flux fluctuations ($`\pm \pi `$) through $`𝐀^s`$, effectively suppressing the holon condensation in (3). Thus in this regime, we expect $`T_v`$ to be reduced to $`T_c`$, as shown in Fig. 3. In particular, since $`a_c`$ also determines the equal time spin-spin correlation length scale, at $`\delta =0.30`$, one has $`a_c`$ $`a`$, which should set an upper doping limit for the bosonic RVB phase, as $`\mathrm{\Delta }^s`$ is associated with nearest-neighbor antiferromagnetic correlations. Both $`T_v`$ and $`T_c`$ are expected to vanish at the point where the bosonic RVB order parameter $`\mathrm{\Delta }^s`$ disappears on the $`T=0`$ axis.

The phase diagram of LSCO has been mapped out carefully, based on high resolution Nernst experiments. For doping concentrations between 0.03 and 0.07, the experimental $`T_v`$ (denoted by $`T_\nu `$ in Ref.) increases very steeply from 0 to $`90`$ K with a slope $`>2,400`$ K, which puts $`(m_ha^2)^1>0.13`$ eV in the present theory. The experimental $`T_v`$ then peaks around $`128`$ $`K`$ at $`\delta =0.11`$ and decreases monotonically at larger doping, a trend similar to the plot shown in Fig. 3.

It is noted that the sharpness of the resonance-like peak at $`E_g`$ is also caused by the holon condensation in the present theory. Thus, $`T_v`$ is intrinsically related to the spin “pseudogap” temperature $`T^{},`$ below which the resonance peak starts to sharpen up in neutron-scattering measurements and the spin-lattice relaxation rate $`1/T_1`$ begins to deviate from the high-T non-Korringa behavior in NMR measurements. In the transport channel, the “pseudogap” seen in resistivity should be also related to the holon condensation. Similarly, the holon condensation plays an essential role in the single-particle channel which may explain the “pseudogap” seen in photoemission spectroscopy measurements. But it is important to point out since no true phase transition takes place and these temperature scales do not necessarily coincide with each other. They represent crossover temperatures in different channels in response to the holon condensation or the forming of the pairing amplitude $`\mathrm{\Delta }_0`$. For instance, the transport $`T^{}`$ is generally in a higher curve, indicating the onset of holon coherence before the holon condensation. If this picture holds true in the cuprates, then one also expects $`T_v`$ to behave like the shaded curve in Fig. 3, which quickly decreases towards $`T_c`$ such that the spontaneous vortex phase shrinks as the doping concentration $`\delta `$ approaches the optimal and over doped regimes.

In summary, we have investigated the normal state below $`T_v,`$ the onset temperature for the amplitude of Cooper pairs, based on the bosonic RVB theory. Such a phase is characterized by spinons that behave as free vortices. $`T_v`$ coincides with the holon condensation, while $`T_c`$ is lower at which the phase coherence is realized due to the binding of spinon-vortices and -antivortices. In the spontaneous vortex phase, the transport properties are quite unique. We showed that free spinon vortices contribute to the Nernst effect, which evolves smoothly into the superconducting phase, consistent with the Nernst experiments. We also argued that the holon condensation is responsible for the pseudogap phenomena and that $`T_v`$ is correlated with various “pseudogap” temperature scales observed experimentally. The bosonic RVB theory of the $`tJ`$ model offers a framework to unify many of the seemingly disparate phenomena observed in the high $`T_c`$ superconductors, and we hope to address these in the future.

###### Acknowledgements.

We thank P. W. Anderson and N. P. Ong for very useful discussions. Z.Y.W. thanks the hospitality of Physics Department at Princeton University where this work was finished, and acknowledges partial support from NSFC Grant 90103021. V. N. M. acknowledges partial support from NSF Grant DMR 98-09483. We thank N. P. Ong and Yayu Wang for sharing their unpublished results with us.

## A Spin rotational invariance

According to the holon effective Hamiltonian (3), the bosonic spinons are perceived by holons as $`\pi `$-flux-tubes, represented by $`𝐀^s`$ defined in (4). Since in the definition (4), the sign of the flux depends on the spin index, one may naturally raise the question as whether the spin rotation may be broken in such a system. In the following we give an explicit proof that the spin rotational symmetry is retained.

For convenience, we will use the original lattice version of $`H_h`$:

$$H_h=t_h\underset{ij}{}\left(e^{i[A_{ij}^e+A_{ij}^s\varphi _{ij}^0]}\right)h_i^{}h_j+H.c.$$
(A1)
where $`A_{ij}^e`$ represents the external electromagnetic field, and the lattice version $`A_{ij}^s`$ is defined by

$$A_{ij}^s=\frac{1}{2}\underset{li,j}{}\text{Im}\mathrm{ln}\left[\frac{z_iz_l}{z_jz_l}\right]\left(\underset{\sigma }{}\sigma n_{l\sigma }^b\right)$$
(A2)
which describes fictitious fluxoids bound to spinons satisfying

$$\underset{c}{}A_{ij}^s=\pm \pi \underset{lc}{}\left[n_l^bn_l^b\right]$$
(A3)
for a closed loop $`c`$. $`\varphi _{ij}^0`$ corresponds to a uniform $`\pi `$ flux per plaquette.

In the phase string formulation, the spin operators are defined as follows:

$$S_i^z=\underset{\sigma }{}\sigma b_{i\sigma }^{}b_{i\sigma },$$
(A4)
$$S_i^+=b_i^{}b_i(1)^ie^{i\mathrm{\Phi }_i^h},$$
(A5)
and $`S_i^{}=(S_i^+)^{}.`$ Here

$$\mathrm{\Phi }_i^h=\underset{li}{}Im\mathrm{ln}\left[z_iz_l\right]n_l^h.$$
(A6)
It is obvious that

$$\left[H_{h,}\text{ }S_l^z\right]=0.$$
(A7)
On the other hand, one finds that the phase $`\mathrm{\Phi }_i^h`$ defined in (A6) plays a crucial role in compensating the extra phase generated from $`A_{ij}^s`$ by a spin flip, which results in

$$\left[H_{h,}\text{ }S_l^\pm \right]=0.$$
(A8)
Also note that if the summation in $`H_h`$ involves the links $`ij`$ which coincide with $`l,`$ such terms have no contribution due to the no-double-occupancy constraint.

So generally one has

$$\left[H_{h,}\text{ }𝐒_i\right]=0$$
(A9)
which means that the holon effective Hamiltonian $`H_h`$ in (A1) or (3) not only always respects the global spin rotational symmetry, but also respects a local spin symmetry. This is consistent with the notion of spin-charge separation in general, and that the spin index is independent of the vorticity of a spinon-vortex composite, discussed in Sec. IIA, in particular. In other words, $`\pi `$-flux-tubes bound to spinons, seen by holons, and the current-vortices bound to spinons only reflect the fact of electron fractionalization, without involving any spin symmetry-breaking. To understand this, one has to realize that there is a peculiar formulation of $`S_i^\pm `$ \[(A5)\] in the phase string formalism, as shown above.


format: text

text:
# Heating of the intracluster medium by quasar outflows

## 1 Introduction

Clusters and groups of galaxies contain a large amount of hot gas, besides galaxies and the gravitationally dominant dark matter. This hot X-ray-emitting gas known as the intracluster medium (ICM) represents a part of the baryonic matter of the universe that is not associated with individual galaxies but remains trapped in the deeper gravitational potential of galaxy clusters. Hierarchical models of structure formation have been very successful in explaining many observed properties of galaxies and galaxy clusters. Nevertheless, some puzzling problems remain open and unexplained. Models of cluster formation in which the intergalactic gas simply falls into the dark matter dominated gravitational potential well (so-called infall models) fail to reproduce all the structural properties of the local cluster population (e.g., Evrard & Henry 1991; Navarro, Frenk & White 1995; Mohr & Evrard 1997; Bryan & Norman 1998). There is certainly some additional physics driving ICM evolution.

Recent X-ray observations have provided evidences for some non-gravitational heating of the diffuse, high density baryons in the potential wells of groups and clusters of galaxies, in addition to the heating during the gravitational collapse. One of the first evidences was in the shape of the $`L_xT`$ relation, which is steeper than the self-similar behaviour $`L_xT^2`$ predicted in the case of gravitational processes only. As early as the emergence of $`ROSAT`$ and $`Einstein`$ data, several authors proposed that the missing element is the existence of a “preheated high-entropy” intergalactic gas prior to a cluster’s collapse (David et al. 1991; Evrard & Henry 1991; Kaiser 1991; White 1991). Later, Ponman et al. (1999), and Lloyd-Davies et al. (2000) found direct evidence of an entropy excess with respect to the level expected from gravitational heating in the centre of groups. Ignoring the constant and logarithms, one can define the “entropy” as $`ST/n^{2/3}`$. The excess entropy, or equivalently, the excess specific energy, flattens the density profile decreasing the X-ray luminosity which is proportional to the square of the density. The effect is stronger in poorer clusters, where the excess energy associated with the excess entropy is comparable to the gravitational binding energy, while rich clusters, where gravity is dominant, are mostly unaffected. This produces a steepening of the $`L_xT`$ relation.

The most popular scenario to successfully explain these thermal properties of the ICM has been the “preheating” scenario. For this scenario, the candidate processes which have been looked into are strong galactic winds driven by supernovae. However Valgeas & Silk (1999) showed that the energy provided by supernovae cannot raise the entropy of intergalactic medium (IGM) up to the level required by current observations. The observed amount of required energy injection depends on the epoch, and have been in the range of 0.4 - 3 KeV per gas particle (Navarro et al. 1995; Cavaliere et al. 1997; Balogh et al.1999; Wu et al. 2000; Lloyd-Davies et al. 2000, Borgani et al. 2001). For example, Wu et al. (2000) showed that galactic winds can impart only $`0.1`$ keV per particle. Moreover, Kravtsov & Yepes (2000) estimated the energy provided by supernovae from the observed metal abundance of ICM and found that the heating only by supernovae driven outflows requires unrealistically high efficiency. On the other hand, quasar outflows may be much more powerful and plausible candidates of the heating (Valageas & Silk 1999). In this paper, we focus on the role of quasar outflows in this regard.

The epoch of the energy input also remains uncertain. Lloyd-Davies et al. (2000) put an upper limit of $`z_{max}7\text{}10`$ on the preheating epoch, from their estimate of excess entropy in groups. For AGNs, there have been no additional constraints like the metal abundance in the case of supernova heating (Kaiser & Alexander 1999). Recently, Yamada & Fujita (2001) have looked into the deformation of the cosmic microwave background (the Sunyaev-Zel’dovich effect) by hot electrons produced at the shocks produced by jets from AGNs. They showed that the observed excess entropy of ICM and $`COBE/FIRAS`$ upper limit for Compton $`y`$-parameter are compatible with each other only when the heating by jets occurred at relatively small redshift ($`z3`$). Thus they questioned the “preheating” scenario as their result suggests that the heating occurred simultaneously with or after cluster formation.

In this paper, we calculate the heat input from quasar outflows inside clusters. We calculate the mechanical work done by various kinds of quasar outflows and the excess energy imparted by them onto the intracluster medium via pdV work. For the statistics of quasars inside clusters, we use the extended Press-Schechter formalism. Finally we calculate the excess energy per particle and tally them with available observations.

In the next section, we discuss the abundance of quasars inside clusters of a given mass. We then discuss the evolution of quasar outflows and the mechanical work done by them in §3. We use these concepts to calculate the heating of the ICM in §4. We then discuss the implications of our results in §5. Throughout the paper we assume a flat universe with a cosmological constant, with $`\mathrm{\Omega }_0=0.35,\mathrm{\Omega }_{\mathrm{\Lambda }0}=0.65,`$ and $`h=0.65`$.

## 2 Quasars inside clusters

For a proper evaluation of the heat input from quasars inside clusters, one first needs to calculate their abundance and its dependence on the quasar mass, cluster mass, and the cluster formation redshift. Observationally, it is still difficult to obtain good statistics of quasars inside clusters. Estimation of the galaxy-QSO correlation function have shown that at low redshifts ($`z0.4`$), quasars typically reside in small to moderate groups of galaxies and not in rich galaxies (e.g., Bahcall & Chokshi 1991; Fisher et al. 1996). It is still uncertain whether or not their optical or radio luminosities depend on their environments. Some studies (e.g., Ellingson, Yee & Green 1991) have shown that radio-loud quasars preferentially reside in richer clusters at higher redshift, although more recent studies (e.g., Wold et al. 2001) do not find any such correlation.

There has been, however, considerable work in relating the observed quasar luminosity function or the radio luminosity function (for radio-loud quasars) with the mass function of galaxies as prescribed by the Press-Schechter (PS) formalism (e.g., Haehnelt & Rees 1993; Haiman & Loeb 1998; Yamada, Sugiyama & Silk 1999).

At $`z=0`$, Yamada, Sugiyama & Silk found that one can reproduce the abundance of the radio sources powered by AGNs (i.e., leaving aside the radio sources powered by starbursts) by assuming that a fraction $`f_r`$ of the halos from PS formalism become radio-loud quasars, where, $`f_r0.01`$ for $`M_h10^{12}`$ M, and $`f_r=0`$ for $`M<10^{12}`$ M. They assume an upper limit on $`M_h`$ of $`10^{14}`$ M. Since it is known that radio loud quasars constitute a fraction $`0.1`$ of the quasar population (Stern et al. 2000), this means that a fraction $`f_q=10f_r`$ of the halos from PS formalism become quasars. In other words, $`f_q0.1`$ for $`10^{12}M_h10^{14}`$ M, and $`f_q=0`$, otherwise. Yamada, Sugiyama & Silk (1999) assumed this fraction to be a constant for all redshifts. In this model, the rate of formation of quasars is given by the derivative of the PS mass function at the relevant mass scale.

This is similar to the model adopted by Haiman & Loeb (1998) and Furlanetto & Loeb (2001; FL01). FL01 showed that at high redshift ($`z4`$) the rate of formation of quasars is a fraction $`f_q0.1`$ of the rate of formation of halos from PS formalism, if a life time of order $`10^7`$ yr is assumed for the quasar. There is, however, a difference, in that they had only a lower limit to $`M_h`$ and no upper limits (see eqns (2) & (3) of Haiman & Loeb 1998). They mention that at low redshifts their formalism does not predict any decline as observed in reality, and for which they consider their model only at high redshift. It is possible that this mismatch is due to the lack of upper limits, since the the differentiation of PS mass function for objects with an upper limit in mass decreases at low redshift, (Haiman, Z. 2001, private communications). In any case, Haiman & Loeb (1998) found that this prescription yields a matching quasar luminosity function that is observed, at redshift $`z2.5`$. We find later that most of the heating of the ICM gas (even for our least massive cluster) occurs at $`z2`$ (Figure 7). We will therefore assume for simplicity that this fraction $`f_q0.1`$ at all redshifts (as in Yamada et al. 1999).

In this paper, we would like to have a conservative estimate of the quasar abundance. Also, we would like to calculate the abundance of quasars in clusters including low mass groups of galaxies. For this reason, we assume the value of $`f_q`$ as above, but use an upper limit of $`10^{13}`$ M. Since the mass function decreases steeply at the higher mass end, this should not change the value of $`f_q`$ substantially. In brief, we assume that,

$$f_q\{\begin{array}{cc}0.1\hfill & \text{if }10^{13}M_h10^{12}M_{}\hfill \\ 0\hfill & \text{if }M_h>10^{13}M_{},M_h<10^{12}M_{},\hfill \end{array}$$
(1)
motivated by the model of Yamada, Sugiyama & Silk (1999), and by the fact that a similar prescription by Haiman & Loeb (1998) recovers the quasar population at high redshift, and relate the rate of formation of quasars with that of halos in the PS formalism.

We are, however, concerned with the statistics of quasars inside clusters. For this one needs to have an extension of the PS mass function which can predict the probability of a given halo becoming a part of bigger object later, or the probability of an object having had a progenitor of a given mass at an earlier epoch. Such extensions of the PS theory have been studied in detail by Bower (1991) and Lacey & Cole (1993), for example.

In the standard PS theory, the mass function, i.e., the fraction of regions with mass in the range $`M,M+dM`$ which have overdensity $`\delta `$ in excess of $`\delta _c`$ (which is the threshold for perturbations becoming non-linear), is given by,

$$f_{PS}=f(M,\delta _c)dM=\frac{1}{\sqrt{2\pi }}\frac{\delta _c}{(\sigma _M^2)^{3/2}}\mathrm{exp}\left[\frac{\delta _c^2}{2\sigma _M^2}\right]\frac{d\sigma _M^2}{dM}dM.$$
(2)
Here, $`\sigma _M`$ is the mass variance of the perturbation at the mass scale $`M`$. The relation between the number density of objects in the mass range $`M,M+dM`$ with $`f(M)dM`$ is,

$$n(M)dM=\frac{\rho _o}{M}f(M)dM,$$
(3)
where $`\rho _o`$ denotes the background mass density.

In the extended PS theory, the fraction of regions of mass $`M`$, contained within a larger scale region of mass $`M^{}`$ and overdensity $`\delta ^{}`$, which are more overdense than $`\delta _c`$, is given by,

$`f(M,\delta _c|M^{},\delta ^{})dM`$ $`={\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle \frac{(\delta _c\delta ^{})}{(\sigma _M^2\sigma _M^{}^2)^{3/2}}}`$ (4)
$`\times \mathrm{exp}\left[{\displaystyle \frac{(\delta _c\delta ^{})^2}{2(\sigma _M^2\sigma _M^{}^2)}}\right]{\displaystyle \frac{d\sigma _M^2}{dM}}dM.`$
This expression recovers the simple PS mass function in the limit $`M^{}\mathrm{}`$ and $`\delta ^{}0`$, relevant for the whole universe.

If we then identify $`M^{}`$ with $`M_{cl}`$, the mass of a cluster, and $`\delta ^{}=\delta _c(z_f)`$, the threshold overdensity of the cluster at its formation epoch $`z_f`$, we can then obtain the mass fraction of $`M_{cl}`$ which have been parts of progenitors of a given mass range ($`M`$ to $`M+dM`$) at a given (earlier) redshift. If we also identify this mass range with that of the quasars as in the standard PS theory ($`10^{12}\text{}10^{13}`$ M), and use the fraction of these halos that become quasars ($`f_q`$; eqn 1), we will obtain the mass fraction of the final cluster which have been quasars at some given earlier epoch.

To obtain the rate of formation of these quasars (inside a future cluster of mass $`M_{cl}`$), we should differentiate the above expression. A simple differentiation will, however, not give the correct result, since there will be a negative contribution from the merging of halos out of this mass range. One would get a negative rate of formation of such quasars at some point if the rate at which they disappear beyond this mass limit is not taken into account in a proper manner. This

Consider the abundance of objects in a given mass range $`M,M+dM`$ at two successive epochs $`z_1`$ and $`z_2`$ ($`z_2<z_1`$). The abundance $`f(M)dM`$ at $`z_2`$ will be given by $`f_2=f_1+FD`$, where $`f_2`$ and $`f_1`$ are the abundances at epochs $`z_2`$ and $`z_1`$, the term $`F`$ denotes the abundance of newly formed objects in this mass range during the epoch $`z_1`$ and $`z_2`$ (from merger of smaller objects), and $`D`$ signifies the abundance of objects that moved out of this mass range as a result of merger (into bigger objects). A simple differentiation of the PS function, involving the difference $`(f_2f_1)`$ will therefore depend on both $`F`$ and $`D`$. For a given range of mass, $`D`$ is very small at a very early epoch, but it increases with time (see, e.g., Fig. 5 of Haiman & Menou 2000), and at a later epoch can become larger than $`F`$. Therefore, at lower redshifts, a simple differentiation can imply a negative rate of change of abundance. If the contribution of $`D`$ is neglected, one would then incorrectly get a negative value of $`F`$.

This problem has been encountered in the case of ordinary PS function by many authors. While studying the rate of mergers in the context of background radiation from starbursts and AGNs, Blain & Longair (1993) noted that a simple differentiation of the PS function leads to a negative rate of formation of objects in a given mass range. In other words, the actual rate of formation of objects is given by,

$$\dot{f}_{form}=\dot{f}_{PS}+\dot{f}_{merger}.$$
(5)
They performed a simulation assuming a simple power spectrum and obtained a fit for the rate of these objects merging to form bigger objects. They found that $`\dot{f}_{merger}`$ can be approximated well (in the Einstein de-Sitter universe) by,

$$\dot{f}_{merger}=\varphi \frac{f_{PS}}{t}\mathrm{exp}\left[(1\alpha )\frac{\delta _c^2}{2\sigma ^2}\right]$$
(6)
where the value of $`\alpha 1.35`$ and $`\varphi 1.3\text{}1.7`$. This problem has also been investigated by Sasaki (1994) and Percival & Miller (1999). With the extension of the PS formalism, one can now calculate this merging rate (see also, Chiu & Ostriker 2000).

We note here that in the case of the merger of a lower mass quasar into a more massive quasar, our implicit assumption is that the central black holes also merge and form a bigger black hole appropriate for the bigger quasar (see below; eqn 13). In other words, we assume that the central black hole mass always traces the halo mass.

The rate at which an object of mass $`M`$ at epoch $`z`$ merges to form a bigger object of mass $`M^{}`$ is given by (Bower 1991; Lacey & Cole 1993),

$`{\displaystyle \frac{d^2p}{dM^{}dz}}(MM^{}|z)dM^{}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi }}}{\displaystyle \frac{\sigma _M^2}{\sigma _M^{}^2(\sigma _M^2\sigma _M^{}^2)^{3/2}}}`$ (7)
$`\times \mathrm{exp}\left[{\displaystyle \frac{\delta _c(z)^2(\sigma _M^2\sigma _M^{}^2)}{2\sigma _M^2\sigma _M^{}^2)}}\right]`$
$`\times |{\displaystyle \frac{d\sigma _M^2}{dM}}||{\displaystyle \frac{d\delta _c(z)}{dz}}|dM.`$
Here, $`p(MM^{}|z)dM^{}`$ is the probability of an object of mass $`M`$ merging to become an object of mass within the range $`M^{},M^{}+dM^{}`$ at redshift $`z`$. The rate of disappearance of objects of a given mass, $`\dot{n}_{merger}`$, should be essentially,

$$\frac{df_{merger}}{dz}(M,z)=f_{PS}(M,z)_{2M}^{\mathrm{}}\frac{d^2p}{dM^{}dz}(MM^{}|z)𝑑M^{},$$
(8)
where the lower limit of the integration is chosen to be such that the merged object is at least twice as massive as the merging object. We show the result of this integration for a sCDM universe (with a COBE normalized power spectrum) in Figure 1, and show the fit of Blain & Longair with $`\varphi =0.9`$ and $`\alpha =1.35`$. We find that the merging rate is fit by a lower value of $`\varphi `$ than they assumed, although the difference is a factor of order unity. It is possible that this difference is due to the specific assumption in the simulation done by Blain & Longair (1993), e.g., in the power spectrum being a simple power law (Blain, A. 2001, private communication), or it can be a result of the lower limit ($`2M`$) chosen by us. At any rate, if we chose the above integral to represent the merger rate then it would be a conservative estimate, since decreasing the lower limit would simple increase the value of the merger rate, and in turn, the formation rate of objects.

We have found that the result in the case of the $`\mathrm{\Lambda }`$CDM universe can be fit by a similar function, with $`\varphi `$ being replaced by $`0.9\frac{d\delta _c(z)}{dz}`$, for $`\alpha =1.35`$ with an accuracy of order $`5\%`$. We, however, do not use these fits in our calculation, and evaluate the integral numerically for our purpose.

To be precise, this rate of disappearance is valid for the objects following the PS mass function, i.e., for objects which are not already parts of bigger objects. Motivated by the extension of the PS formalism, we here posit that the rate of disappearance of objects inside a bigger object also has the same form, with $`f_{PS}(M,\delta _c)`$ in eqn (8) being replaced by $`f(M,\delta _c|M^{},\delta ^{})`$. There is admittedly no way of verifying the truth of this ansatz at present, since this would involve more extensions of the PS theory than that is available now. It will also involve comparing the merger rates inside and outside of clusters. It, however, leads to a conservative estimate for the formation rate of quasar in a cluster. As the work of Bower (1991) has shown, growth of perturbations inside a cluster is enhanced compared to in the field. This means that the merging rate of objects of a given mass inside a cluster should be larger than that in the field. Here, by assuming a comparable merging probability (the factor that multiplies the abundance of objects, given by $`f_{PS}`$), we are in a way underestimating the rate of disappearance ($`\dot{f}_{merger}`$), and in turn, the rate of formation ($`\dot{f}_{form}`$) of quasars in a cluster. The final result of total heat input from our formalism should, therefore, be a conservative estimate.

We show the results of adding the rate of disappearance in Figure 2. The dashed lines, which show the term $`(df_{PS}/dz)`$ (or equivalenlty, $`\frac{df_{PS}}{dt}`$), become negative at lower redshift, suggesting the need for the addition of the rate of disappearance of objects. The dotted lines show the result of adding this rate, using the Blain & Longair (1993) fit (with $`\varphi =0.9`$), and the solid lines show the result using the integral in equation 8. The upper panels of the figure show the case of objects (of mass $`10^{13}`$ M) in the field, and the bottom panels show the case for these objects inside a cluster of $`M_{cl}=10^{15}`$ M ($`z_f=0`$). For these curves, we have used the integration in equation 8, and the extension of the PS mass function for the abundance of objects inside the cluster, as explained above. The left and right panels show the cases for sCDM and $`\mathrm{\Lambda }`$-CDM universe. As the bottom panels show, the addition of the disappearance rate does not suffice to make $`[(df_{PS}/dz)]`$ (inside clusters) positive at low redshift. This is suggestive of the enhanced growth of perturbation, and the need for a larger rate of merger inside clusters, as previous authors have noted.

We also show the case for $`M_{cl}=10^{14}`$ M by the long dashed line in the lower right panel of Figure 2. Comparison with the solid line (for $`M_{cl}=10^{15}`$ M) shows that the formation rate of these galaxies inside a lower massive cluster is larger. This follows simply from the extension of the PS formalism (from the dependence on the term $`(\sigma _M^2\sigma _M^{}^2)`$). It is interesting to note that this is consistent with the suggestion from observation that quasars are preferentially located in groups of galaxies instead of rich clusters (e.g., Bahcall & Chokshi 1991; Fisher et al. 1996).

If our formalism is used without any correction, this will lead to subtraction of energy input in the final result. We circumvent this problem by putting $`[(df_{PS}/dz)]=0`$ when this term turns negative. This will, therefore, provide a lower limit to the total energy input from quasar outflows in a cluster.

We can finally write down the rate of formation of quasars in a given mass range $`M,M+dM`$ inside a cluster of a given mass, $`M_{cl}`$, in the form of the rate of increase of the fraction of mass of $`M_{cl}`$ which is in quasars at an epoch $`z`$, as (remembering the equation 3),

$`{\displaystyle \frac{df_{q,cl}(M,z|M_{cl},z_f)}{dz}}`$ $`dM=f_q{\displaystyle \frac{df(M,\delta _c(z)|M_{cl},\delta _c(zf))}{dz}}dM`$ (9)
$`+f_qf(M,\delta _c(z)|M_{cl},\delta _c(zf))dM`$
$`\times {\displaystyle _{2M}^{M_{cl}}}{\displaystyle \frac{d^2p}{dM^{}dz}}(MM^{}|z)dM^{},`$
with the condition that $`\frac{dn_q}{dz}=0`$, for $`z<z_n`$, where $`\frac{dn_q}{dz}|_{z<z_n}>0`$. The integral on the right hand side is evaluated using eqn(7). Here we have also changed the upper limit of the integration to $`M_{cl}`$. The integrand is a rapidly decreasing function of $`M^{}`$ and the value of the integral depends mostly on the lower limit.

Here, the threshold density contrast in a cosmological constant dominated universe is given by a fit given by Kitayama & Suto (1996),

$$\delta _c(z)=1.68[g(z=0)/g(z)][1+0.0123\mathrm{log}\mathrm{\Omega }_m(z)].$$
(10)
In our calculations, we have used a fit for $`g(z)`$ from Carroll, Press & Turner (1992),

$`g(\mathrm{\Omega }_m(z),\mathrm{\Omega }_\mathrm{\Lambda }(z))`$
$`{\displaystyle \frac{5\mathrm{\Omega }_m(z)}{2[\mathrm{\Omega }_m(z)^{4/7}\mathrm{\Omega }_\mathrm{\Lambda }(z)+(1+\frac{\mathrm{\Omega }_m(z)}{2})(1+\frac{\mathrm{\Omega }_\mathrm{\Lambda }(z)}{70})]}},`$ (11)
where (Lahav et al. 1991),

$`\mathrm{\Omega }_m(z)`$ $`=`$ $`\mathrm{\Omega }_{m0}(1+z)^3/[\mathrm{\Omega }_{m0}(1+z)^3+\mathrm{\Omega }_{\mathrm{\Lambda }0}],`$
$`\mathrm{\Omega }_\mathrm{\Lambda }(z)`$ $`=`$ $`\mathrm{\Omega }_{\mathrm{\Lambda }0}/[\mathrm{\Omega }_{m0}(1+z)^3+\mathrm{\Omega }_{\mathrm{\Lambda }0}].`$ (12)
We emphasize here that the above formalism leads to a conservative estimate of abundance of quasars in a cluster (and, therefore, the final heat input), because (a) we ignore the increased pace of growth of perturbation and the merging rate inside a cluster, and (b) the lower limit of the integration could in reality be smaller than $`2M`$, which is probably the reason the rate of formation still turns negative at low redshifts even after the addition of merger term.

To summarise the work in this section, we have used the existing ideas for relating the quasar formation rate to the Press-Schechter mass function, to estimate the rate of formation of quasars inside clusters (as a function of cluster mass and formation redshift), utilizing the extensions of PS formalism. First, we use eqn (1) to relate the PS mass function to quasar abundance. The standard PS mass function (eqn 2) is then replaced by its extension (eqn 4). Furthermore, we add the contribution due to merger (into larger objects) (eqns 5 & 8), again using the extensions of PS mass function. We finally have the rate of formation of quasars inside clusters as given by eqn (9).

We use this formalism to calculate the total energy input from outflows from quasars in a cluster, or, a lower limit to it. We next discuss the energy input from individual outflows, which we will combine with our calculation of formation rate of quasars for the final result.

## 3 Work done by quasar outflows

In this section we calculate the energy input from quasar outflows into the ambient medium. We model the outflows as they evolve in the ambient medium and calculate the pdV work done by the outflows. To begin with, we discuss different kinds of outflows that we consider and the characteristics of the hosts of quasars.

### 3.1 Quasar outflows

We consider two major types of quasar outflows. For radio-loud quasars (RLQ), the outflow is in the form of a tightly collimated jet, which deposits energetic particles into a cocoon which expands against the surrounding medium. These outflows are characterised by the kinetic luminosity of the jet, $`L_k`$. According to Willot et al. (1999), this is correlated with the bolometric luminosity $`L_{bol}`$ of the quasar, and that $`0.05L_k/L_{bol}1.0`$. We follow FL01 in arguing that since $`L_{bol}10L_B`$ (Elvis et al. (1994)), the rest-frame B-band luminosity, $`L_kL_B`$.

Radio-loud quasars, however, constitute only about $`10\%`$ of the total population of quasars (Stern et al. 2000). We therefore define a factor $`f_o`$ for the fraction of quasars with outflows, and define $`f_o0.1`$ for our RLQ model. The fraction $`f_o`$ here denotes the number of radio relics/lobes per halo, since a radio loud quasar may have several outbursts of radio activity. This fraction is therefore a very conservative estimate since it is obtained from observed radio luminosity function and does not take into account the existence of radio relics in clusters.

Another important kind of outflows are encountered in broad absorption line (BAL) quasars. The absorption troughs are thought to be due to absorbing clouds flowing out of the quasars with velocities up to $`0.1c`$. Although they are encountered in about $`10\%`$ of the quasars, it is believed that all quasars have such outflows (all the time) and the covering fraction of the BAL outflows is about $`10\%`$ (Weymann et al. 1991; Weymann 1997). Some authors also believe that BAL outflows have a limited lifetime (especially the low ionization BALs) and that they have a large covering fraction in the early phase of a quasar (Voit et al. 1993). For our calculation, we use a fraction $`f_o1`$ for the BAL outflows. We discuss the effect of the uncertainty in these factors on the final result in §5.

We model the BAL outflows as having a kinetic luminosity $`L_k`$. Following FL01, if $`N_H`$ is the column density of the absorbing gas, $`f_c`$ the covering fraction, and $`R_{BAL}`$ is the size of the absorption system, then $`L_k`$ is related to the outflow velocity $`v_{BAL}`$ as $`L_k2\pi f_cN_Hm_pR_{BAL}v_{BAL}^3`$. The observed range of these parameters are as follows: $`v_{BAL}0.1c`$, $`f_c0.1`$ (Weymann 1997; but see above), $`R_{BAL}1\text{}500`$pc and $`N_H10^{22}\text{}10^{23}`$ cm<sup>-2</sup> (Krolik 1999; Gallagher et al. 1999). For these values, the magnitude of $`L_k`$ is close to that of $`L_B`$. FL01 also argued that for BAL winds $`L_k0.1\text{}100L_B`$, and finally assumed $`L_kL_B`$. Since this estimate depends crucially on a number of uncertain parameters (for example, the fact that the absorption column density in optical measurement is much smaller than the above mentioned X-ray column density), it may not really be a conservative estimate, but it does provide a simple scaling which we hope is not too unreasonable. The estimate $`L_KL_B0.1L_{Edd}`$ is probably not a conservative estimate, but an upper limit, in that for a covering fraction of $`10\%`$, the mechanical luminosity of the wind could not be larger than $`0.1`$ of the Eddington rate. Keeping all these uncertainties in mind, we assume that $`L_kL_B`$ for BAL outflows.

We then need to connect $`L_B`$ of a given quasar with the properties of its halo. Firstly, as Haiman & Loeb (1998) have shown, the mass of the black hole at the centre is related to $`L_B`$, as,

$$M_{BH}=\frac{1}{0.093}\left(\frac{L_B}{1.4\times 10^{38}\mathrm{erg}\mathrm{s}^1}\right)M_{}$$
(13)
assuming that the quasar radiates at the Eddington luminosity. The factor of $`0.093`$ reflects the fraction of the Eddington luminosity radiated in the B-band, taken from the median quasar spectrum of Elvis et al. (1993). Statistically speaking, we therefore assume that a fraction $`f_q`$ (eqn(1)of all black holes radiate at the Eddington rate for a life time of $`10^7`$ yr, while the rest does not radiate at all.

Secondly, the correlation between the central black hole mass and the total baryonic mass of the galaxy (Magorrian et al. 1998; Gebhardt et al. 2000) gives $`M_{BH}4\times 10^4M_h`$, where $`M_h`$ is the total mass of the galaxy, using a value of $`M_{BH}/M_{baryonic}2\text{}3\times 10^3`$, and $`M_{baryonic}/M_h\mathrm{\Omega }_b/\mathrm{\Omega }_00.2`$.

As far as the collimation is concerned, the geometry of BAL outflows is still uncertain, whereas the radio jets are well collimated. Since some models do suggest a modest collimation even in BAL outflows, with a covering fraction of $`f_c0.1`$ (Weymann 1997), we use the idea of collimated outflows for outflows from both radio-loud and BAL quasars.

We next discuss the evolution of the outflows from radio-loud quasars and calculate the fraction of its total kinetic luminosity that it deposits into the surrounding medium in the form of pdV work. For concreteness, we will assume that the corresponding fraction for BAL outflows has similar values, and keep in mind the uncertainty in the geometry and energetics of BAL outflows.

### 3.2 Evolution of outflows

The standard scenario for outflows from radio loud quasars involves a ‘cocoon’ surrounding the core and the jet, and consisting of a shocked ambient medium and shocked jet material (Scheuer 1974; Blandford and Rees 1974). Begelman & Cioffi (1989) constructed a simple model of the evolution of a cocoon in which the cocoon is overpressured against the ICM. In their model, the expansion along the jet axis is determined by the balance of the thrust of the jet and the ram pressure, whereas the thermal pressure of the cocoon drives along the direction perpendicular to the jet axis. Results of numerical simulations agree with this scenario (Loken et al. 1992; Cioffi & Blondin 1992).

Here we adopt the model of the evolution of cocoons following the approach of Bicknell et al. (1997), which is based on the Begalman & Cioffi (1989) model but includes the pdV work done by the cocoon, in order to find the fraction of total energy lost by the quasar to the ICM through mechanical work (pdV work). Bicknell et al. (1997) derived this fraction ($`f_{pdV}`$) to be $`f_{pdV}=3/8`$, for a homogeneous ambient medium (their equation (2.13)).

We will, however, calculate this fraction from numerical solution of the equations governing the evolution of the cocoon, for the following reasons. Firstly, the derivation mentioned above implicitly assumes that the mean pressure averaged over the hotspot region is equal to the mean lobe pressure (in the language of Bicknell et al. (1997), this means $`\zeta 1`$). In fact, their equation (2.13) shows that for constant $`\zeta `$, one has in general, $`P_c\frac{dV_c}{dt}=(1+2\zeta )L_j/8\zeta `$, which recovers the fraction $`3/8`$ for $`\zeta =1`$. In reality, however, this ratio does not remain a constant in time. Secondly, this derivation is valid only during the period when the jet is active. Even after the jet switches off, the cocoon, however, continues to evolve as a result of its overpressure until it reaches an equilibrium pressure with the ambient medium (see also Nath 1995). The cocoon, therefore, continues to do pdV work even after the jet switches off, and the inclusion of this process will lead to an upward revision of the fraction of total energy that is lost in pdV work. Besides, there seems to be some confusion in the literature regarding the fraction. For example, Inoue & Sasaki (2001) have recently adopted a fraction $`f_{pdV}=1/4`$ in their calculation of energy input into the surrounding gas.

We therefore calculate this fraction by numerically solving the equations of cocoon evolution.

We consider two collimated steady jets advancing into the ambient ICM. The thermalized jet matter and the shock-compressed ICM matter form a cocoon around the jets and the cocoon expands with shocks advancing in directions both parallel and perpendicular to the jet axis. After this stage of evolution, when the jet turns off after a lifetime of $`t_{life}3\times 10^7`$ years (Kaiser 2000), cocoons still retain high pressure. They cool radiatively and expand due to its overpressure till it reaches a pressure equilibrium with the ambient medium. Thus the relevant equations are:

$`{\displaystyle \frac{dr_h}{dt}}`$ $`=`$ $`\left({\displaystyle \frac{L_j}{A_h\rho _a\beta c}}\right)^{1/2},t<t_{life}`$ (14)
$`=`$ $`\left({\displaystyle \frac{P_c}{\rho _a}}\right)^{1/2},t>t_{life}`$
$`{\displaystyle \frac{dr_c}{dt}}`$ $`=`$ $`\left({\displaystyle \frac{P_c}{\rho _a}}\right)^{1/2},`$ (15)
$`{\displaystyle \frac{dE_c}{dt}}`$ $`=`$ $`L_jP_c{\displaystyle \frac{dV_c}{dt}}`$ (16)
where $`L_j`$ is the jet luminosity, $`\rho _a`$ is the density of the ambient medium and $`\beta c`$ is the velocity of the jet material. As the jet is highly relativistic, $`\beta 1`$. The averaged hotspot area $`A_h30`$ kpc<sup>2</sup> (Bicknell et al. 1997) is assumed to be larger than the radius of the jet, according to the ‘dentist’s drill’ model of the jet (Scheuer 1982). Here $`r_h`$ is the length of the jet or the distance of the hotspot from the centre of the galaxy, $`r_c`$ is the half-width of the cocoon at the centre and $`V_c`$ is the volume of the cocoon given by $`V_c=ϵ_v(2\pi r_c^2)r_h`$ where $`ϵ_v`$ is the geometrical factor depending on the shape of the cocoon. For our calculations we have taken the shape of the cocoon as biconical and so $`ϵ_v1/3`$. Finally $`P_c`$ is the pressure inside the cocoon given by $`P_c=(E_c/V_c)(\gamma 1)`$ where $`\gamma =4/3`$ and $`E_c`$ is the total energy inside the cocoon given by $`E_c=L_jt`$ till $`t<t_{life}`$ and $`E_c=L_jt_{life}`$ afterwards, where $`t_{life}`$ is the lifetime of the jet.

This is admittedly a simplified model of the evolution of the cocoon. In reality, after the jet switches off, one expects Rayleigh-Taylor and Kelvin-Helmholtz instabilities to distort the cocoon, giving rise to ‘buoyant’ plumes. This phase of the evolution of the cocoon, and its effect on the ambient medium, has been recently addressed by various authors (Gull & Northover 1973; Churazov et al. 2000; Brüggen & Kaiser 2000), mainly with the help of numerical simulations. It is possible that this phase adds substantial heating to the intracluster medium, but short of doing a numerical simulation it is difficult to assess its importance. We have also neglected the loss of energy through radiation, as modeling such losses would involve detail knowledge of various parameters (e.g., electron energy spectrum and the possibilities of re-acceleration of electrons). With the uncertainties involved in modeling these process, it seems reasonable to adopt the above simplified picture as a pointer and keep the uncertainties in mind while discussing the final result. In light of this discussion, we will also calculate the final result with a value of $`f_{pdV}=3/8`$ as in Bicknell et al. (1997), which we will adopt as a conservative lower limit.

We numerically calculate the volume of the cocoon as it grows into the ICM and the pressure inside the cocoon at each step, and add up the pdV work to get the final amount of energy lost in this mode. The evolution of the cocoon is continued until the pressure inside the cocoon becomes equal to the ambient pressure, $`n_ak_BT_a`$ ($`n_a=\rho _a/\mu m_p`$), where $`T_a`$ is the temperature of the ICM. The fraction $`f_{pdV}`$ is calculated by taking the ratio of the total energy lost through mechanical work to the total energy (that is $`L_jt_{life}`$). (We found that for the relevant values of the ambient medium parameters, the time scale to reach pressure equilibrium is always larger than $`t_{life}`$.) The dependence of $`f_{pdV}`$ on the ambient density, is shown in Figure 2 for $`T_a=10^6,10^7`$ and $`10^8`$ K, and for $`L_k=10^{46}`$ (solid lines) and $`10^{47}`$ erg/s (dashed lines).

The plot shows that the fraction $`f_{pdV}`$ is a function of temperature of the cluster and also the density of the ambient medium. The general trend is that the fraction reduces at higher temperatures and higher densities. This is because of the fact that the cocoon reaches pressure equilibrium with the ambient medium sooner for a higher pressure environment (higher $`T_a`$ and $`n_a`$), and the total pdV work ends up being smaller. The plot also shows that the fraction $`f_{pdV}`$ depends weakly on the jet luminosity for lower temperatures ($`T_a10^7`$ K), whereas there is a bit of a difference for $`T_a10^8`$ K.

The fraction $`f_{pdV}`$ calculated above is somewhat larger than that has been used in the literature, for ambient medium with low pressure. The difference is mainly the result of our inclusion of cocoon evolution even after the jet has switched off. Incidentally, Inoue & Sasaki (2001), while using a value of $`f_{pdV}=1/4`$, discussed the possibility that this fraction could be larger in reality, because of its continued evolution after the switching off of the jet (their §3.2).

In our calculation for the total pdV work done by BAL and RLQ outflows, we will use the values of $`f_{pdV}`$ obtained above. As mentioned earlier, the energetics and geometry of BAL outflows are not clear at present. For concreteness, we have worked out the case of RLQ outflows in detail, and we will use the same values of $`f_{pdV}`$ for BAL outflows as well.

## 4 Heating of the ICM

Equipped with the knowledge of the rate of formation of quasars in clusters (equation 9) and the fraction of total energy which is deposited as pdV work by the outflows from them (§3.2), we are now in a position to calculate the total amount of non-gravitational energy provided by quasar outflows in a cluster. If we denote the gas fraction of the total cluster mass by $`f_{gas}`$, then the total number of gas particles is $`M_{cl}f_{gas}/m_p`$, where $`m_p`$ is the proton mass. It is not yet clear if the gas fraction has a universal value for clusters of all masses and at all redshifts. A few authors (see, e.g., Schindler (1999)) have found no correlation of the gas fraction with the cluster mass, whereas others (e.g., Ettori & Fabian (1999)) find some correlation (with low mass clusters possessing lower gas fraction). This correlation has been attributed to the excess energy deposition (since low mass clusters are more liable to lose gas from heating than massive clusters) (see, e.g., Bialek et al. (2000)). Here, we are, however, trying to calculate the magnitude of this very excess energy. It would not be appropriate to include this correlation it a priori in our calculation. (We show later that including these correlations only increases our estimate of excess energy input.) We have, therefore, used a value of $`f_{gas}=0.1`$ for our calculation. The total energy per (gas) particle deposited into the ICM of a cluster of mass $`M_{cl}`$ is then given by,

$`ϵ_{pdV}=`$ $`{\displaystyle \frac{m_p}{M_{cl}f_{gas}}}{\displaystyle _{z_m}^0}{\displaystyle _{M_l}^{M_u}}{\displaystyle \frac{df_{q,cl}(M,z|M_{cl},z_f)}{dMdz}}`$ (17)
$`{\displaystyle \frac{M_{cl}}{M}}dMdzf_o[L_kt_{life}f_{pdV}(n_a,T_a)],`$
where $`\left[\frac{d^2n_q(M,z|M_{cl},z_f)}{dMdz}dMdz\right]`$ is calculated using equation 9. The factor $`f_o1`$ for the BAL outflows, and $`f_o0.1`$ for outflows from RLQs (§3). The redshift $`z_m`$ is the maximum redshift of heat input. We later show (Figure 7) that the heat input is negligible for $`z5`$. The density and temperature of the ICM of a cluster of a given mass ($`M_{cl}`$) and formation redshift ($`z_f`$) is calculated using (Eke et al. 1998),

$$T_a=1.65\times 10^7(1+z)\left[\frac{M_{cl}}{10^{15}h^1M_{}}\right]^{2/3}\left[\frac{\mathrm{\Omega }_0\mathrm{\Delta }(\mathrm{\Omega }_0,z)}{\mathrm{\Omega }(z)}\right]^{1/3},$$
(18)
and,

$$n_a=\frac{M_{cl}f_{gas}}{m_p(4/3)\pi r_{vir}^3},$$
(19)
where

$$\mathrm{\Delta }(z)=18\pi ^2+82x39x^2$$
(20)
and $`x=\mathrm{\Omega }(z)1`$ (Bryan & Norman 1998), where we use equation 12 to compute $`\mathrm{\Omega }(z)`$. For $`r_{vir}`$, we use,

$$r_{vir}=\left(\frac{3M_{vir}}{4\pi \mathrm{\Delta }(z)\rho _{crit}(z)}\right)^{1/3},$$
(21)
where we have used $`M_{vir}=M_{cl}`$, and $`\rho _{crit}`$ is the critical density of the universe. The densities used in the following calculations range between $`10^4\text{}10^6`$ cm<sup>-3</sup>.

The integral in equation (17) is evaluated using $`M_l=10^{12}`$ M and $`M_u=10^{13}`$ M, for $`z_f=0,0.5,1`$ for different values of $`M_{cl}`$. We present the results for the total non-gravitational energy input per particle as a function of cluster mass (or, equivalently, gas temperature) in Figure 4 (for $`z_f=0`$). The solid curve shows the heat input calculated using $`f_{pdV}`$ from §3.2. The dotted line shows the case for a constant $`f_{pdv}=3/8`$ (Bicknell et al. 1997). We show the results for different $`z_f`$ in Figure 5 (against $`T`$) and Figure 6 (against cluster mass).

We also show in Figure 7 the rate of deposition of energy ($`\frac{dϵ_{pdV}}{dz}`$) as functions of the redshift for three clusters of masses $`M=2\times 10^{13},10^{14}`$ and $`10^{15}`$ M, all for $`z_f=0`$. All the curves drop to zero at low redshift because of the condition $`dn_{q,cl}/dt=0`$ ( in eqn 9). In reality the contribution to the heating should be small but non-zero, and will increase the estimate of excess energy.

## 5 Discussion

It has been estimated that the amount of excess energy required to explain the observation is of order $`0.5\text{}3`$ keV per particle (Navarro et al. 1995; Cavaliere et al. 1997; Balogh et al.1999; Wu et al. 2000). Recently, however, Lloyd-Davies et al. (2000) have shown from observations of groups of clusters that an excess energy of $`0.44\pm 0.3`$ keV per particle suffices to explain the excess entropy in groups. They showed that this can explain the entropy floor for galaxy groups with temperature $`T4`$ keV. Borgani et al. (2001) have also shown with the aid of numerical simulations that excess energy of order $`1`$ keV per particle reproduces the observations.

The solid and dashed curves in Figures 4 & 5 show that the excess energy from pdV work done by quasar outflows fall in this required range. It is seen that the excess energy per particle is larger for clusters or groups with lower temperature. This is due to two factors: (a) the number of quasars per unit mass is larger for smaller clusters, and (b) the fraction of total energy in outflows that is lost in pdV work is larger for them. We show the results in the case of a constant $`f_{pdV}=3/8`$ (as in Bicknell et al. (1997)) with dotted lines. It is interesting to note that even in this case the excess energy is in the required range ($`0.44\pm 0.3`$ keV per particle), especially for groups with low temperatures, as advocated by Lloyd-Davies et al. (2000). Incidentally, this is larger than the estimate of excess energy from galactic winds ($`0.1`$ keV per particle, Wu et al. (2000)).

We found that our results for the excess energy (solid line in Figure 4) can be approximated by a fit of type (in keV per particle),

$`ϵ_{pdV}`$ $`0.258\left({\displaystyle \frac{T}{10\mathrm{keV}}}\right)^{0.193}+0.033\left({\displaystyle \frac{T}{2\mathrm{keV}}}\right)^{1.2},`$ (22)
$`0.50T8.0\mathrm{keV}`$
We also found that the results for the excess energy taking $`f_{pdV}=3/8`$ (Bicknell et al. (1997)) (dotted line in Figure 4) can be approximated by a fit of type (in keV per particle),

$$ϵ_{pdV}0.17+0.045\left(\frac{T}{2\mathrm{keV}}\right)^{0.95},0.50T8.0\mathrm{keV}$$
(23)
We compare our results with the data from Lloyd-Davies et al. (200) in Figure 8. We show the predictions of our calculation as the thin solid curve (corresponding to the solid curve in Figure 4), where the data points have been taken from that of Lloyd-Davies et al. (2000) (their Figure 9). We also show the result of the calculations for $`f_{pdV}=3/8`$ by the thick solid line. The data points refer to the binding energy of a (constant) central fraction ($`0.004`$) of the virial mass of groups and clusters. The dashed line shows the case for self-similar models, of type $`ET`$, derived from the data points for rich clusters. The dotted line shows their fit (with a constant excess energy of $`0.44`$ keV per particle) along with a formal 1 $`\sigma `$ confidence interval shown by the shaded region. The figure shows that our predictions are consistent with the data available at present. The thick line (corresponding to $`f_{pdV}=3/8`$) falls close to the fit provided by Lloyd-Davies et al. (2000), whereas the thin line (using $`f_{pdV}`$ from Figure 3) somewhat overestimates the heat input at the low mass end. The thick line can be viewed as a conservative estimate of the heat input, since it uses $`f_{pdV}=3/8`$. The thin line, however, provides an estimate of the heat input if $`f_{pdV}`$ is much larger than $`3/8`$. We should remind ourselves here that we have calculated $`f_{pdV}`$ for radio galaxies and used the same values for BAL outflows. If a more accurate estimate of $`f_{pdV}`$ for BAL outflows is worked out in the future, the resulting heat input into the ICM could then be scaled accordingly using Figure 8.

We would like to emphasize here again that the our calculation provides a conservative estimate of the excess energy, for reasons outlined in §2. Moreover, we have used a constant density and temperature in time for the ICM gas (for clusters with a given $`z_f`$), which is not very realistic. In reality, the density at higher redshift will be smaller, and the inclusion of a density dependent $`f_{pdV}`$ will only increase the estimate of excess energy (since this fraction increases with decreasing density).

The curves for clusters forming at different epochs show excess energies to decrease somewhat for clusters with higher formation redshift.

The curves of Figures 4 & 5 assume $`f_o1`$, which is relevant for BAL outflows. The excess energy from RLQ outflows will be one tenth of these curves, showing the difficulty of using radio galaxies as the only source of non-gravitational heating, if conservative estimates for their kinetic luminosities are used. Recently, Inoue & Sasaki (2001) have used the radio luminosity functions of Willot et al. (2001) and Ledlow & Owen (1996) to determine the abundance of radio galaxies in clusters, and finally to estimate the total pdV work done by the cocoons of these radio galaxies. They estimated an excess energy of order $`1`$ keV per particle for rich clusters like the Coma cluster, and also for poor groups, assuming that their ratio of radio galaxies per unit cluster mass is universal. From our calculation, we find an excess energy from only radio-loud quasars that is an order of magnitude smaller than their estimate. It is possible that the assumptions leading to the estimate of $`L_k`$ are at the source of this difference (see, e.g., the discussion on the uncertainty in the factor $`f_j`$ in their §3.2).

The evolution of the rate of energy input (Figure 7) shows that the heat input before $`z5`$ is almost negligible. Most of the heating occurs in the range $`z1\text{}4`$. It is also clearly seen that the ICM of poor groups is heated at lower redshifts, compared to the gas in massive clusters. This follows from the simple consideration that the evolution of objects in a given mass range (here, that of the quasars) typically occurs earlier in more massive clusters. Our results, therefore, suggests that the ICM in clusters were “preheated”, before the major mergers took place in them, whereas, the ICM in groups of galaxies were heated at epochs similar to that of their formation. It is interesting to note that the energy input epoch lies below the upper limit from recent observations of Lloyd-Davies et al. (2000). It is also consistent with the limit on the range of redshifts as shown by Yamada & Fujita (2001) given the level of uncertainty.

Recently, accelerated particles from shocks as a result of the formation of clusters have been hypothesized to be the source for the diffuse gamma ray background (Loeb & Waxman 2000). This gamma ray production will be suppressed, however, if the gas in the clusters were heated substantially at earlier epochs (Totani & Inoue 2001). Our results show that the ICM in massive clusters were preheated (and earlier than the ICM in groups), and the suppression of the gamma ray production will, therefore, be an important effect, if confirmed.

Finally, we discuss the uncertainties involved in our calculation. Apart from the uncertainties in cosmological parameters, the major uncertainties lies in the factors $`f_q`$ (§2) (connecting the abundances of quasars with PS mass function) , $`f_o`$ (the fraction of quasars with outflows) ,$`f_{pdV}`$ and $`f_{gas}`$. Among these, the most uncertain factor is $`f_o`$, which we have assumed to be of the order of unity for BAL outflows. The uncertainty in this factor will be reflected in the uncertainty of the final heat input (with a direct proportionality; see eqn 17). The uncertainty in $`f_{pdV}`$ has been already discussed earlier, and we found that even if $`f_{pdV}`$ is as low as $`3/8`$ for all cases, the final excess energy is certainly larger than that from supernovae driven winds, and is still within the required range of excess energy, especially for loose groups. As far as the uncertainty in $`f_{gas}`$ is concerned, we have also done our calculation with a varying $`f_{gas}`$, e.g., of the type,

$$f_{gas}=0.15(1+z)^{0.5}(M_{cl}/10^{15}h^1)^{0.1},$$
(24)
as has been advocated by Ettori & Fabian (1999), and we have found that the excess energy is approximately doubled in this case. It is, however, not clear if this correlation is a result of the excess energy, and, so, it would not be appropriate to attach much significance to this result. We have also varied the lower limit in our estimate of $`f_q`$ (eq. 1) and found that changing the lower limit from $`10^{12}`$ M to $`10^{11}`$ M increases the final heat input by only $`10`$ %. This is because of the fact that the increase in the number of quasars is compensated by the decrease in their mechanical luminosity. Lastly, we have already discussed in detail the uncertainty in the net formation rate of quasars in clusters, and as explained in §2, our approach here has been very conservative, and the final results should be regarded as conservative estimates in this regard.

## 6 Summary

We have calculated the excess energy deposited by quasar outflows in clusters in order to explain the observations of excess entropy in groups and clusters of galaxies. We summarise our findings below:

(1) We have used the extended Press-Sechter formalism to derive a formation rate of quasars inside clusters and groups, as a function of the cluster/group mass and its formation redshift.

(2) We have calculated the fraction of the kinetic luminosity of outflows (RLQ and BAL outflows) that is deposited onto the ambient medium, as a function of the density and temperature of the ambient medium. For outflows from radio-loud quasars, we have included the evolution of the cocoon after the jet turns off.

(3) The final excess energy from the mechanical work done by quasar outflows is found to be of order $`0.18\text{}0.85`$ keV per particle, and is consistent with the data available at present. The excess energy in this scenario comes mainly from BAL outflows, with radio galaxies supplying about a tenth of the total. Keeping in mind the uncertainties in the estimate of energetics and abundances of radio and BAL outflows, we conclude that both radio galaxies and BAL outflows are promising candidates for heating the ICM. We found that this excess energy increases with decreasing mass of the cluster/group. This prediction could be tested with better data in the near future. The excess energy does not depend strongly on the formation redshift.

(4) The epoch of heating is found to be in the range $`z1\text{}4`$, where this epoch is at lower redshifts for low mass clusters.

The authors would like to thank Nahum Arav, Andrew Blain, Sergio Colafrancesco, Zoltan Haiman, Subhbrata Majumdar and Trevor Ponman for valuable discussions. We have also benefitted from the detail comments of the anonymous referee. We are grateful to Ed Lloyd-Davies for supplying the data points for binding energy. We would also like to thank Rekhesh Mohan for help with plotting packages.