text:
# Critical temperature and condensate fraction of the trapped interacting Bose gas with finite-size effects

## I Introduction

The experimental realization of Bose-Einstein condensation (BEC) in the dilute alkali-metal atoms  and more recently in the atomic hydrogen  has stimulated a new interest in the theoretical study of the inhomogeneous Bose gas. Thermodynamic properties such as critical temperature, condensate fraction have been investigated by several authors for the trapped Bose gas recently .

In the presence of an anisotropic harmonic potential of the form $`V_{ext}(\stackrel{}{r})=m(\omega _x^2x^2+\omega _y^2y^2+\omega _z^2z^2)/2`$ the noninteracting model gives the critical temperature $`T_c^0=\frac{\mathrm{}\omega _{ho}}{k_B}\left(\frac{N}{\zeta (3)}\right)^{1/3}`$ , where $`\omega _{ho}=\left(\omega _x\omega _y\omega _z\right)^{1/3}`$ is the geometrical average of the oscillator frequencies. $`\zeta (n)`$ is Riemann $`\zeta `$ function. In the large-N limit, the condensate fraction is $`\xi =\frac{N_0}{N}=`$ $`1\frac{\zeta (3)}{N}\left(\frac{k_BT}{\mathrm{}\omega _{ho}}\right)^3`$. Finite-size effects  and interatomic interaction will give corrections to the thermodynamic properties of the system. The correction to the critical temperature due to finite-size effects has shown to be $`\delta T_c^0/T_c^00.73\overline{\omega }N^{1/3}/\omega _{ho}`$, where $`\overline{\omega }=(\omega _x+\omega _y+\omega _z)/3`$ is the mean frequency.

Although the atom clouds realized in the present experiments are very dilute, the effects due to interatomic interaction are important at low temperature. Researches show that atom-atom interaction will give leading corrections when $`N`$ is larger than $`10^5`$ for the alkalis. In fact the question of how two-body forces affect the thermodynamic properties of these systems have been the object of several theoretical investigations . Using a local density approximation Giorgini, Pitaevskii, and Stringari  obtained the shift of the critical temperature due to interatomic interaction: $`\delta T_c^0/T_c^01.33aN^{1/6}/a_{ho}`$, where $`a_{ho}=\sqrt{\mathrm{}/m\omega _{ho}}`$ is the harmonic oscillator length. Recently Monte Carlo simulation  is also used to investigate the thermodynamic properties of the trapped interacting Bose gases.

In this work, we investigate the critical temperature and the condensate fraction of the system when both atom-atom interaction and finite-size effects are taken into account. Canonical ensemble is used to obtain the analytical high-order correction due to interatomic interaction. Especially, finite-size effects are obtained simultaneously. To obtain the shift of the critical temperature, we give the analytic description of the condensation fraction near the critical temperature.

The paper is planed as follows. In Sec.II we outline the canonical ensemble. In Sec.III we give the condensate fraction of the system near the critical temperature. The shift of the critical temperature agrees with the usual results . In Sec.IV we obtain the analytical high-order correction to the condensate fraction due to atom-atom interaction. The lowest-order correction agrees with the well-established results . In Sec.V the correction to the condensate fraction beyond mean field theory is given.

## II Partition Function of the Trapped Interacting Bose Gases

Canonical ensemble and saddle-point method have been used to investigate the thermodynamic properties of the interacting uniform Bose gases . In this work canonical ensemble is used to discuss the trapped interacting Bose gases. The partition function of $`N`$ trapped interacting bosons is given by

$$Q(N)=\mathrm{\Sigma }_n\mathrm{exp}\left(\beta E_n\right)=\mathrm{\Sigma }_{\{n_x,n_y,n_z\}}\mathrm{exp}\left[\beta \left\{\mathrm{\Sigma }_{n_x,n_y,n_z}N_{n_xn_yn_z}\epsilon _{n_xn_yn_z}+E_{int}\right\}\right],$$
(1)
where $`N_{n_xn_yn_z}`$ and $`\epsilon _{n_xn_yn_z}`$ are the occupation numbers and energy level of the state $`\{n_x,n_y,n_z\}`$ respectively. $`E_{int}`$ is the interaction energy of the system.

From (1)

$$Q(N)=\mathrm{\Sigma }_{N_0=0}^N\left\{\mathrm{exp}\left[\beta \left(E_0+E_{int}\right)\right]Q_0\left(NN_0\right)\right\},$$
(2)
where $`Q_0\left(NN_0\right)=\mathrm{\Sigma }_{\{n_x,n_y,n_z\}}^{}\mathrm{exp}\left[\beta \left(\mathrm{\Sigma }_{n_x,n_y,n_z}N_{n_xn_yn_z}\epsilon _{n_xn_yn_z}\right)\right]`$ stands for the partition function of a fictitious system of $`NN_0`$ trapped noninteracting atoms. $`E_0`$ is the energy of the condensate. For convenience, we have separated out the ground state $`n_x=n_y=n_z=0`$ from the states $`\{n_x,n_y,n_z\}0`$, so that we first carry out the primed summation over all $`\{n_x,n_y,n_z\}`$ ($`\{n_x,n_y,n_z\}0`$) with a fixed value $`(NN_0)`$ of the partial sum $`\mathrm{\Sigma }_{\left\{n_xn_yn_z\right\}}^{}`$ and then carry out a summation over all possible values of $`N_0`$, ie. $`N_0=0`$ to $`N_0=N`$. Assume $`A_0\left(NN_0\right)`$ is the free energy of the fictitious system. $`A_0(NN_0)=\frac{1}{\beta }\mathrm{ln}Q_0(NN_0)`$. From (2)

$$Q(N)=\mathrm{\Sigma }_{N_0=0}^N\mathrm{exp}\left[\beta \left(E_0+E_{int}\right)\beta A_0(NN_0)\right].$$
(3)
The sum $`\mathrm{\Sigma }_{N_0=0}^N`$ in (3) maybe replaced by the largest term in the sum, for the error omitted in doing so will be statistically negligible. With this approximation we do not investigate the fluctuations of the system. However, this approximation is reasonable because we can give the occupation number of the energy level $`\epsilon _{n_xn_yn_z}`$, which agrees with the widely used mean occupation number in the frame of grand-canonical ensemble. Neglecting all terms but the largest in (3), the number of bosons in the condensate can be obtained as:

$$\beta \frac{}{N_0}\left(E_0+E_{int}\right)\beta \frac{}{N_0}A_0(NN_0)=0.$$
(4)
The calculations of the free energy $`A_0(NN_0)`$ of the fictitious noninteracting Bose gas is nontrivial because there is a requirement that the number of the particles is $`NN_0`$ in the summation of the partition function. Using the saddle-point method of integration developed by Darwin and Fowler  it is straightforward to find that $`\beta \frac{}{N_0}A_0(NN_0)=\mathrm{ln}z_0`$ in which $`z_0`$ is the fugacity of the fictitious Bose gas. In addition, the fugacity $`z_0`$ is determined by

$$NN_0=\mathrm{\Sigma }_{n_x,n_y,n_z}^{}\frac{1}{\mathrm{exp}\left[\beta \epsilon _{n_xn_yn_z}\right]z_0^11}.$$
(5)
We can easily understand (5) in terms of Bose-Einstein distribution of the trapped noninteracting Bose gas. Using the relation $`\beta \frac{}{N_0}A_0(NN_0)=\mathrm{ln}z_0`$, (4) becomes

$$\beta \frac{}{N_0}\left(E_0+E_{int}\right)+lnz_0=0.$$
(6)
Equations (5) and (6) will be used to discuss the corrections due to atom-atom interaction and finite-size effects. Once we know the interaction energy of the system, it is easy to obtain the correction due to atom-atom interaction. In addition, finite-size effects are separated off in (5). Omitting interactions between atoms, we can obtain $`\mathrm{ln}z=\beta ϵ_{000}`$ from (6). From (5) the number of the thermal atoms $`N_T`$ is given by

$$N_T=NN_0=\mathrm{\Sigma }_{n_x,n_y,n_z}^{}\frac{1}{\mathrm{exp}\left[\beta (\epsilon _{n_xn_yn_z}ϵ_{000})\right]1}.$$
(7)
This is the exact conclusion in the frame of grand-canonical ensemble for the trapped noninteracting Bose gas. This shows the equivalence between the grand-canonical ensemble and canonical ensemble in the calculations of the condensate fraction. In addition, it shows that neglecting all terms but the largest in (3) is reasonable. In the absence of interaction, from (7), the condensate fraction is given by

$$\xi =1(\frac{T}{T_c^0})^3\frac{3\overline{\omega }\zeta \left(2\right)}{2\omega _{ho}\zeta \left(3\right)^{2/3}}(\frac{T}{T_c^0})^2N^{1/3}.$$
(8)
The third term in (8) is the finite-size correction to the ideal Bose gas in the thermodynamic limit. In short, (5) and (6) account for both interatomic interaction and finite-size effects. We will use (5) and (6) to calculate the critical temperature and condensate fraction of the system.

## III critical temperature and condensate fraction near $`T_c`$

### A interaction energy of the system near the critical temperature

The parameter expressing the importance of the interatomic interaction compared to the kinetic energy is $`\frac{E_{int}}{E_{kin}}\frac{N_0|a|}{a_{ho}}`$, where $`a`$ is the scattering length between bosons. Near the critical temperature $`\frac{N_0|a|}{a_{ho}}<<1`$. This means that $`E_{kin}>>E_{int}`$. In this case we can use the method of pseudopotentials developed by Huang and Yang  to calculate the interaction energy of the system.

In the method of pseudopotentials, the actual Hamiltonian of the system is replaced by an effective Hamiltonian, such that the ground state and the low-lying energy levels of the system are given equally well by the new Hamiltonian. With hard-sphere approximation, the boundary conditions between atoms are replaced by the pseudopotential operator $`g\delta ^{(3)}(\stackrel{}{r})\frac{}{r}r`$ in the effective Hamiltonian, where $`g=\frac{4\pi a\mathrm{}^2}{m}`$. The effective Hamiltonian of the system may be taken to be

$$\widehat{H}=\frac{\mathrm{}^2}{2m}\left(_1^2+\mathrm{}+_N^2\right)+\mathrm{\Sigma }_{i=1}^NV_{ext}(\stackrel{}{r}_i)+\mathrm{\Sigma }_{i<j}\omega _{ij},$$
(9)
where $`\omega _{ij}=g\mathrm{\Sigma }_{i<j}\delta ^{(3)}\left(\stackrel{}{r}_i\stackrel{}{r}_j\right)\frac{}{r_{ij}}r_{ij}`$. For a trapped dilute Bose gas near the critical temperature the pseudopotentials can be regarded as perturbation terms. The energy levels to the first order in $`a`$ may be obtained through the usual perturbation theory.

In the general case of a system containing an arbitrary number $`N`$ of bosons, the unperturbed normalized wave function of the system is given by

$$\mathrm{\Psi }_N=\left(\frac{N_{\alpha _1}!N_{\alpha _2}!\mathrm{}}{N!}\right)^{1/2}\mathrm{\Sigma }_P[\varphi _{\alpha _1}\varphi _{\alpha _2}\mathrm{}\varphi _{\alpha _N}],$$
(10)
where the sum is taken over all permutations of the different suffixes $`\alpha _1`$, $`\alpha _2`$, $`\mathrm{}`$, $`\alpha _N`$ and the numbers $`N_{\alpha _i}`$ show how many of these suffixes have the same value $`\alpha _i`$ (with $`\mathrm{\Sigma }N_{\alpha _i}=N`$). From the first-order perturbation theory the interaction energy $`E_{int}`$ of the system is given by

$$E_{int}=<\mathrm{\Psi }_N|\mathrm{\Sigma }_{i<j}\omega _{ij}|\mathrm{\Psi }_N>.$$
(11)
Thus, the operators $`\left(/r_{ij}\right)r_{ij}`$ will be operating on a set of functions which are well behaved for all values of $`r_{ij}`$; accordingly, these operators may be replaced by the unit operators. Basing on the investigation of the various permutations  of the single particle states the interaction energy takes the form,

$$E_{int}=\frac{g(N_0^2N_0)}{2}\varphi _0^4(\stackrel{}{r})d^3\stackrel{}{r}+2gN_0\varphi _0^2(\stackrel{}{r})n_T\left(\stackrel{}{r}\right)d^3\stackrel{}{r}+gn_T^2(\stackrel{}{r})d^3\stackrel{}{r}.$$
(12)
To obtain (12) we used the fact that $`n_T(\stackrel{}{r})=\mathrm{\Sigma }_{l0}N_l\varphi _l^2(\stackrel{}{r})`$ is the density distribution of the thermal atoms. Near the critical temperature, the density distribution of the condensate is $`n_0\left(\stackrel{}{r}\right)=N_0\varphi _0^2\left(\stackrel{}{r}\right)=N_0\left(\frac{m\omega _{ho}}{\pi \mathrm{}}\right)^{3/2}exp[\frac{m}{\mathrm{}}\left(\omega _xx^2+\omega _yy^2+\omega _zz^2\right)]`$. The density distribution of the thermal atoms can be obtained using Bose-Einstein distribution and the semiclassical approximation of the energy level of a single particle. Near the critical temperature $`\frac{k_BT}{\mathrm{}\omega _{ho}}N^{1/3}`$. In the available traps $`N`$ ranges from a few thousand to several millions, thus $`\mathrm{}\omega _{ho}<<k_BT_c^0`$, ie. the energy level of the thermal atoms can be approximated as continuous. In addition, this means that the semiclassical approximation for the normal gas is expected to work well on a wide range of temperatures . In terms of Bose-Einstein distribution $`n_T\left(\stackrel{}{r}\right)=d^3\stackrel{}{p}\left(2\pi \mathrm{}\right)^3\left[\mathrm{exp}\left(\beta \epsilon (\stackrel{}{p},\stackrel{}{r})\right)1\right]^1`$, where $`\epsilon (\stackrel{}{p},\stackrel{}{r})=\frac{\stackrel{}{p}^2}{2m}+V_{ext}\left(\stackrel{}{r}\right)`$ is the semiclassical energy in phase space and $`\beta =1/k_BT`$. Density distribution of the thermal atoms in the coordinate space would be $`n_T(\stackrel{}{r})=\lambda _T^3g_{3/2}\left(\mathrm{exp}[V_{ext}(\stackrel{}{r})/(k_BT)]\right)`$, where $`\lambda _T=\left[2\pi \mathrm{}^2/\left(mk_BT\right)\right]^{1/2}`$ is the thermal wavelength. $`g_{3/2}\left(z\right)`$ belongs to the class of functions $`g_\alpha (z)=\mathrm{\Sigma }_{n=1}^{\mathrm{}}z^n/n^\alpha `$. By integrating the semiclassical approximation $`n_T(\stackrel{}{r})`$ over space one obtain the condensate fraction $`\xi =1\left(\frac{T}{T_c^0}\right)^3`$, which does not account for finite-size effects. However, in this paper the semiclassical approximation $`n_T(\stackrel{}{r})`$ is used to calculate the interaction energy of the system, ie. used to calculate the high-order correction due to interatomic interaction. Thus the adoption of the semiclassical approximation will omit only higher-order modification, which is of the order of the multiplication of the corrections due to the interatomic interaction and finite-size effects. Because the finite-size effects have been separated off in (5), we can give finite-size effects although semiclassical approximation is used to calculate the correction due to interatomic interaction.

Using the geometrical average of the oscillator frequencies $`\omega _{ho}`$, $`n_T\left(r\right)=\lambda _T^3g_{3/2}\left(\mathrm{exp}\left[m\omega _{ho}^2r^2/2k_BT\right]\right)`$. Introducing the width of the normal gas $`R_T=\sqrt{2k_BT/m\omega _{ho}^2}`$ and defining the rescaled variable through $`r=R_T\overline{r}`$, we can obtain $`n_T\left(\overline{r}\right)=\lambda _T^3g_{3/2}\left(\mathrm{exp}\left[\overline{r}^2\right]\right)`$. Near the critical temperature $`k_BT>>\mathrm{}\omega _{ho}`$, the classical Boltzmann distribution $`n_{cl}\left(\overline{r}\right)\lambda _T^3\mathrm{exp}\left[\overline{r}^2\right]`$ gives a zero-order approximation for the density distribution of the thermal cloud. Thus we take the trial function of the form

$$n_T^G\left(\overline{r}\right)=\eta _1\lambda _T^3\mathrm{exp}\left[\eta _2\overline{r}^2\right]$$
(13)
where $`\eta _1`$ and $`\eta _2`$ are two dimensionless variational parameters. Another reason for the adoption of the Gaussian trial function lies in the fact that there is a factor $`\mathrm{exp}\left[\overline{r}^2\right]`$ in $`n_T\left(\overline{r}\right)`$. $`\eta _1`$ and $`\eta _2`$ are determined by the requirement that $`|n_T^G\left(\stackrel{}{r}\right)n_T\left(\stackrel{}{r}\right)|d^3\stackrel{}{r}`$ obtain the minimum. Results of numerical calculations show that $`\eta _1=2.587`$, $`\eta _2=1.146`$ and $`|n_T^G\left(\stackrel{}{r}\right)n_T\left(\stackrel{}{r}\right)|d^3\stackrel{}{r}/n_T\left(\stackrel{}{r}\right)d^3\stackrel{}{r}=0.0017`$. So the Gaussian distribution $`n_T^G(\stackrel{}{r})`$ agrees very well with $`n_T(\stackrel{}{r})`$. We should note that the Gaussian distribution $`n_T^G\left(\stackrel{}{r}\right)`$ is different from the classical Boltzmann distribution because of the variational parameters $`\eta _1`$ and $`\eta _2`$.

Combining (12) and (13) we have

$$E_{int}=\frac{g}{2\left(\sqrt{2\pi }a_{ho}\right)^3}N_0^2+\left(2\frac{1}{2^{3/2}}\right)\eta _1g\lambda _T^3N_0+\frac{\eta _1g\lambda _T^3}{2^{3/2}}N,$$
(14)
where we have used the fact that $`a_{ho}>>\lambda _T`$ and $`k_BT>>\mathrm{}\overline{\omega }`$ near the critical temperature.

### B shift of the critical temperature and condensate fraction near $`T_c`$

We introduce a scaling parameter which accounts for the role of the two-body repulsive interaction. This parameter $`\theta `$ is fixed by the ratio between $`gn_T\left(\stackrel{}{r}=0,T_c^0\right)`$ and the critical temperature $`k_BT_c^0`$. $`\theta =\frac{gn_T\left(\stackrel{}{r}=0,T_c^0\right)}{k_BT_c^0}=\frac{\eta _1\sqrt{2/\pi }}{\zeta \left(3\right)^{1/6}}\frac{a}{a_{ho}}N^{1/6}`$. The scaling parameter $`\theta `$ can also be written in the form $`\theta =0.65\eta ^{5/2}`$, where $`\eta =\mu \left(T=0\right)/k_BT_c^0=\frac{\zeta \left(3\right)^{1/2}15^{2/5}}{2}\left(N^{1/6}\frac{a}{a_{ho}}\right)^{2/5}`$ is also an important scaling parameter accounting for the role of interatomic interaction . Obviously $`\theta `$ stands for high-order correction due to interatomic interaction, compared to the parameter $`\eta `$. From (5), (6), and (14) we can obtain the equation on the condensate fraction.

$$1\xi =\frac{t^3}{\zeta \left(3\right)}g_3\left(z_0\right)+\frac{3\zeta \left(2\right)}{2\zeta \left(3\right)^{2/3}}\frac{\overline{\omega }}{\omega _{ho}}t^2N^{1/3},$$
(15)
$$\frac{\zeta \left(3\right)^{1/2}\theta N^{1/2}}{\eta _1t}\xi \left(2\frac{1}{2^{3/2}}\right)\theta t^{1/2}+\mathrm{ln}z_0=0,$$
(16)
where we have used the reduced temperature $`t=T/T_c^0`$. From (15) and (16) the condensate fraction is given by

$$\xi =\frac{1t^34.51\frac{a}{a_{ho}}N^{1/6}t^{7/2}\frac{3\zeta \left(2\right)}{2\zeta \left(3\right)^{2/3}}\frac{\overline{\omega }}{\omega _{ho}}t^2N^{1/3}}{1+1.16\frac{a}{a_{ho}}N^{2/3}t^2}.$$
(17)
The third term in the numerator of (17) represents the correction due to the interaction between atoms, while the last term in the numerator accounts for the correction due to the finite-size effects. By setting $`\xi =0`$ in (17) one can obtain the shift of the critical temperature,

$$\frac{\delta T_c^0}{T_c^0}=\frac{\delta T_{int}}{T_c^0}+\frac{\delta T_{finite}}{T_c^0}=1.50\frac{a}{a_{ho}}N^{1/6}\frac{\zeta \left(2\right)}{2\zeta \left(3\right)^{2/3}}\frac{\overline{\omega }}{\omega _{ho}}N^{1/3}.$$
(18)
The first term in (18) is the shift due to interatomic interaction. It agrees with the results based on the local density approximation . Because of the denominator in (17), close to $`T_c^0`$ the condensate fraction will increase slowly with decreasing the temperature. Especially, it shows that there is only high-order correction of the scaling parameter $`\eta `$ due to interatomic interaction. We should note that (17) holds only when $`\frac{N_0a}{a_{ho}}<<1`$. The second contribution in (18) gives exactly the usual results due to the finite-size effects.

## IV condensate fraction of the system below $`T_c`$

### A lowest-order correction due to interatomic interaction and finite-size effects

Below the critical temperature $`\frac{N_0a}{a_{ho}}>>1`$, $`\mathrm{𝑖𝑒}.`$ $`E_{int}>>E_{kin}`$. In this case we can use the well-known Thomas-Fermi approximation. With the Thomas-Fermi approximation the density profile of the condensate is $`n_0\left(\stackrel{}{r}\right)=\frac{\mu V_{ext}\left(\stackrel{}{r}\right)}{g}\theta \left(\mu V_{ext}\left(\stackrel{}{r}\right)\right)`$. The normalization condition on $`n_0(\stackrel{}{r})`$ provides the relation between chemical potential and number of particles in the condensate: $`\mu =\frac{\mathrm{}\omega _{ho}}{2}\left(\frac{15N_0a}{a_{ho}}\right)^{2/5}`$. From Gross-Pitaevskii (GP) equation  and the well-known Virial theorem the energy of the condensate turns out to be $`E_0=\left(5/7\right)\mu N_0`$ . Omitting atom-atom interactions in the normal gas and the interaction between the condensate and normal gas, from (5) and (6)

$$N_0=NN_T=N\mathrm{\Sigma }_{n_x,n_y,n_z}^{}\frac{1}{\mathrm{exp}\left[\beta \left(\epsilon _{n_xn_yn_z}\mu \right)\right]1}.$$
(19)
From (19) it is easy to obtain the following result

$$\xi =1t^3\frac{\zeta \left(2\right)}{\zeta \left(3\right)}\frac{\mu }{k_BT}\frac{3\overline{\omega }\zeta \left(2\right)}{2\omega _{ho}\zeta \left(3\right)^{2/3}}t^2N^{1/3},$$
(20)
where the last term in (20) is the correction due to finite-size effects . Using the relation $`\mu /k_BT=\eta \xi ^{2/5}/t`$, the third term in (20) gives exactly the usual lowest-order modification due to atom-atom interactions .

$$\xi =1t^3\frac{\zeta \left(2\right)}{\zeta \left(3\right)}\eta \xi ^{2/5}t^2\frac{3\overline{\omega }\zeta \left(2\right)}{2\omega _{ho}\zeta \left(3\right)^{2/3}}t^2N^{1/3}.$$
(21)
This proves the equivalence between canonical ensemble and grand-canonical ensemble, even in the presence of interatomic interactions. We will give high-order modification in the following.

### B high-order modification due to interatomic interactions

Equation (21) is obtained by omitting the interaction between thermal atoms and interaction between thermal atoms and the condensate. When these two sorts of interactions are considered

$$E_0+E_{int}=\frac{5}{7}\mu N_0+2gn_0\left(\stackrel{}{r}\right)n_T\left(\stackrel{}{r}\right)d^3\stackrel{}{r}+gn_T^2\left(\stackrel{}{r}\right)d^3\stackrel{}{r}.$$
(22)
The calculations of (22) are straightforward. The results is given by

$$E_0+E_{int}=\frac{5\mu N_0}{7}+\frac{16\pi \eta _1\mu }{15}\frac{R_\mu ^3}{\lambda _T^3}\frac{16\pi \eta _1\eta _2}{35}\frac{R_\mu ^3}{\lambda _T^3}\frac{\mu ^2}{k_BT}+\frac{\eta _1g\lambda _T^3}{2^{3/2}}\left(NN_0\right),$$
(23)
where $`R_\mu =\sqrt{2\mu /m\omega _{ho}^2}`$ is the radius of the condensate. $`R_\mu ^3/\lambda _T^3=\alpha _{R\lambda }N^{1/2}t^{3/2}N_0^{3/5}`$, where $`\alpha _{R\lambda }=0.294\left(a/a_{ho}\right)^{3/5}`$. From (23)

$$\beta \frac{}{N_0}\left(E_0+E_{int}\right)=\left(10.89\eta ^{5/2}t^{1/2}\right)\frac{\mu }{k_BT}1.01\eta ^{5/2}t^{1/2}.$$
(24)
From (5), (6), and (24) we can obtain the condensate fraction below the critical temperature,

$$\xi =1t^3\frac{\zeta \left(2\right)}{\zeta \left(3\right)}t^3\left[\left(10.89\eta ^{5/2}t^{1/2}\right)\frac{\eta \xi ^{2/5}}{t}+1.01\eta ^{5/2}t^{1/2}\right]\frac{3\overline{\omega }\zeta \left(2\right)}{2\omega _{ho}\zeta \left(3\right)^{2/3}}t^2N^{1/3}.$$
(25)
In (25) the terms comprises $`\eta ^{5/2}`$ represents high-order correction to the condensate fraction due to interatomic interaction. Omitting the terms comprises $`\eta ^{5/2}`$, we can obtain the lowest-order correction (21). Essentially, (21) and (25) are transcendental equations on the condensate fraction. We can easily give the numerical solution of the condensate fraction. In Fig.1 and Fig.2 we used the experimental parameter by Ensher et al. , where the cloud consists of $`4\times 10^4`$ atoms at the transition and $`a/a_{ho}=5.4\times 10^3`$. According to Fig.2 the high-order correction (25) agrees well with the condensate fraction (17) in the region $`\frac{N_0a}{a_{ho}}<<1`$. From (17) and (25), close to $`T_c^0`$ the condensate fraction increase quite slowly with decreasing the temperature. However, the lowest-order correction (21) is rather different from the dashed line near the critical temperature. The difference between (17) and (21) lies in the fact that there is only high-order correction due to interatomic interaction near the critical temperature, while there is only lowest-order correction in (21). Fig.2 demonstrates the validity of (25) over the whole range of temperature below $`T_c`$.

In the preceding calculations we neglects the role of the interaction between the condensate and the thermal atoms on the density distributions of the condensate and the thermal atoms. In the frame of Hartree-Fork model , the densities of the condensate and the thermal component are given as follows.

$$n_0\left(\stackrel{}{r}\right)=\frac{\mu V_{ext}\left(\stackrel{}{r}\right)2gn_T\left(\stackrel{}{r}\right)}{g}\theta \left(\mu V_{ext}\left(\stackrel{}{r}\right)2gn_T\left(\stackrel{}{r}\right)\right),$$
(26)
$$n_T\left(\stackrel{}{r}\right)=\lambda _T^3g_{3/2}\left[e^{\left(v_{ext}\left(\stackrel{}{r}\right)+2g\left(n_0\left(\stackrel{}{r}\right)+n_T\left(\stackrel{}{r}\right)\right)\mu \right)/k_BT}\right].$$
(27)
Using (5), (6), (22) and (26), (27) and the usual iterative procedure we can obtain the numerical conclusion of the condensate fraction. The numerical conclusion is illustrated in Fig.3.

## V condensate fraction beyond mean-field theory

In the preceding calculations of the condensate fraction, GP equation is used to obtain the energy $`E_0`$ of the condensate. GP equation is expected to be valid if the system is dilute, ie. $`n|a|^3<<1`$. Obviously, correction to GP equation would modify the condensate fraction of the system. The first correction to the mean-field approximation have been investigated by Timmermans, Tommasini, and Huang  and by Braaten and Nieto . For large $`N_0`$, using the local density approximation, the density distribution of the condensate is given by

$$n_0\left(\stackrel{}{r}\right)=\frac{\mu V_{ext}\left(\stackrel{}{r}\right)}{g}\frac{4m^{3/2}}{3\pi ^2\mathrm{}^2}\left[\mu V_{ext}\left(\stackrel{}{r}\right)\right]^{3/2},$$
(28)
with $`\mu `$ given by

$$\mu =\frac{\mathrm{}\omega _{ho}}{2}\left(\frac{15N_0a}{a_{ho}}\right)^{2/5}\left(1+\sqrt{\pi a^3n\left(0\right)}\right).$$
(29)
The parameter $`a^3n\left(0\right)`$ can be directly expressed in terms of the relevant parameters of the system. $`a^3n\left(0\right)=\frac{15^{2/5}}{8\pi }\left(N_0^{1/6}\frac{a}{a_{ho}}\right)^{12/5}=0.0078\eta ^6\xi ^{2/5}`$. Using the relation $`E_0=5N_0\mu /7`$ and (5), (6) we can obtain the correction to the condensate fraction below $`T_c`$

$$\xi =1t^3$$

$$\frac{\zeta \left(2\right)}{\zeta \left(3\right)}t^3\left[\left(10.89\eta ^{5/2}t^{1/2}+0.16\eta ^3\xi ^{1/5}\right)\frac{\eta \xi ^{2/5}}{t}+1.01\eta ^{5/2}t^{1/2}\right]\frac{3\overline{\omega }\zeta \left(2\right)}{2\omega _{ho}\zeta \left(3\right)^{2/3}}t^2N^{1/3}.$$
(30)
The term $`0.16\eta ^3\xi ^{1/5}`$ is the correction due to the modification of the GP equation. It gives higher-order correction of the scaling parameter $`\eta `$, compared to the high-order correction due to interatomic interactions.

Another correction to the condensate fraction beyond mean-field approximation is given by the quantum depletion of the condensate. There are atoms which do not occupy the condensate at zero temperature because of quantum depletion effects. We can use local density approximation to write the density of atoms out of the condensate. Timmermans, Tommasini, and Huang  gives $`n_{out}\left(\stackrel{}{r}\right)=\left(8/3\right)\left[n\left(\stackrel{}{r}\right)a^3/\pi \right]^{1/2}`$. Integration of $`n_{out}\left(\stackrel{}{r}\right)`$ yields the result

$$\xi _{out}=n_{out}\left(\stackrel{}{r}\right)d^3\stackrel{}{r}=\frac{5\sqrt{\pi }}{8}\sqrt{a^3n\left(0\right)}.$$
(31)
Combining (30) and (31)

$$\xi =1t^3\frac{\zeta \left(2\right)}{\zeta \left(3\right)}t^3\left[\left(10.89\eta ^{5/2}t^{1/2}+0.16\eta ^3\xi ^{1/5}\right)\frac{\eta \xi ^{2/5}}{t}+1.01\eta ^{5/2}t^{1/2}\right]$$

$$0.11\eta ^3\xi ^{1/5}\frac{3\overline{\omega }\zeta \left(2\right)}{2\omega _{ho}\zeta \left(3\right)^{2/3}}t^2N^{1/3}.$$
(32)
The term $`0.11\eta ^3\xi ^{1/5}`$ is the correction due to quantum depletion.

Compared to the finite-size effects, the corrections to the mean-field approximation have a different dependence on the parameter $`N`$ and $`a/a_{ho}`$. The corrections beyond mean field theory become larger than finite-size effects when $`N`$ is larger than about $`10^6`$. In Fig.4 we give numerical conclusion of (32). The corrections beyond mean field theory become important when the temperature is much lower than the critical temperature.

## VI conclusion

To conclude, the canonical ensemble is used to obtain the equations on the condensate fraction (5) and (6), which account for interaction effects. Finite-size effects can be obtained simultaneously because these effects are comprised in (5). (5) and (6) are used to investigate the condensate fraction near and below the critical temperature. From the condensate fraction near the critical temperature, we obtain the shift of the critical temperature due to atom-atom interaction and finite-size effects simultaneously. In addition, the analytical high-order correction due to interatomic interaction is obtained in this work.

## Acknowledgments

This work was supported by the Science Foundation of Zhijiang College, Zhejiang University of Technology. This work was also supported by the National Natural Science Foundation of China (19975019).


format: text

text:
# Essential and inessential features of Hawking radiation

## 1 Introduction

Hawking radiation from black holes is a semiclassical quantum effect that has now been with us for some 27 years , and whose theoretical importance is difficult to exaggerate. Over the decades, the Hawking effect has accreted a quite considerable mythology. Perhaps the two most pernicious myths attached to this effect are:

* Hawking radiation has something to do with gravity”, and,
* Hawking radiation automatically implies Bekenstein entropy”.

These myths were engendered by two historical accidents: (1) the Hawking effect was first encountered within the context of general relativity, and (2) the fact that it was discovered shortly after the notion of Bekenstein entropy (geometric entropy) had been formulated .

Though the Hawking effect was at first partly motivated by the need for a consistent thermodynamic interpretation for the notion of Bekenstein entropy, it was rapidly appreciated that the Hawking effect is much more primitive and fundamental  — in particular the Hawking effect continues to make sense even in situations where geometric entropy and even gravity itself are simply not relevant . This observation underlies much of the current interest in “analog models of/for general relativity” ; there is now a realistic possibility for experimental detection of the Hawking effect in condensed-matter analog systems using current or planned technology , with “effective metrics” and “black holes” that have nothing to do with gravity itself . In view of this situation, in this current article I will attempt to isolate an irreducible minimum of physical assumptions needed for the Hawking effect to arise.

Since the literature is vast, I will not be able to do justice to all known derivations of the effect — some key derivations are those due to Hawking himself , the Hartle–Hawking approach using analytic continuation of the propagator across the event horizon of an eternal black hole , and the Gibbons–Hawking approach using Euclideanization (Wick rotation) . Extremely useful general surveys are provided in  and .

A particularly relevant discussion is the early work of Damour and Ruffini , who emphasise the behaviour of the “outgoing” modes as one crosses the horizon. Though their presentation is given in terms of the Kerr–Newman geometry it is easy to verify that the specifics of the geometry enter *only* in the evaluation of the surface gravity of the future Killing horizon, and that the route from surface gravity to Hawking radiation does not depend on either the Einstein equations or even the underlying physical mechanism leading to the existence of the metric.

In a slightly different vein, the discussion of Parker  particularly emphasises the relationship with particle production from a dynamical vacuum state, while Gerlach , Grove , Hu , and Brout and Parentani  emphasise in varying degree the near-universal role of the exponential stretching associated with many types of horizon.

More recently the contributions of Massar and Parentani , Parikh and Wilczek  Padmanabhan et al , and Schutzhold  should be noted. They emphasise, in slightly different forms, the analyticity properties of the modes and what is effectively an imaginary contribution to the action localized at the horizon, an observation that can be traced back to the work of Damour and Ruffini .

It is also important to realise that the Hawking radiation effect is independent of whatever cutoff you introduce to the high-frequency physics — this is one of the theoretical reasons for interest in analog models, because for acoustic black holes you have an explicit model for the high frequency cutoff in terms of atomic physics . I will not have anything specific to say about this particular issue in the present paper.

Based on the results extracted from 27 years of research, as verified and made more explicit by the recent interest in analog models, it is clear that the basic physics we should be aiming for is:

* Hawking radiation is kinematics”,
* Bekenstein entropy is geometrodynamics”.

That is: Hawking radiation is a purely kinematic effect that depends only on the existence of a Lorentzian metric (with no particular prejudice as to how this metric arises) and some sort of horizon; Hawking radiation does not depend on the validity of the Einstein equations (as may most quickly be verified by looking at Hawking’s original derivation  and verifying that the Einstein equations are nowhere used nor needed). In contrast, Bekenstein’s geometric entropy associated with the area of the event horizon is an intrinsically geometrodynamic effect, wrapped up with the validity of the Einstein equations: Entropy equals one quarter the area (plus perturbative corrections) if and only if the Einstein equations are valid (plus perturbative corrections) . The connection with the Einstein equations arises because when integrating the first law to evaluate the Bekenstein entropy you need to use the relationship between total mass of the black hole and its surface gravity (and hence Hawking temperature), and it is in this relationship between surface gravity and mass-energy that the Einstein equations enter.

I shall also distinguish the notion of “apparent horizon” from that of the “event horizon” (absolute horizon) and demonstrate that the existence of a locally definable apparent horizon is quite sufficient for obtaining the Hawking effect. (Remember that to define the event horizon you need to know the entire history of the spacetime out to the infinite future; you should be a little alarmed if the question of whether or not a black hole is radiating *now* depends on what it is doing in the infinite future.)

The general theme of the analysis will be to do as much as possible with the eikonal approximation (even WKB is mild overkill), and look for generic features in the modes at/near the apparent horizon. Specifically, we will look for a Boltzmann factor. Note that a good derivation should not depend on either grey-body factors or a past horizon. Grey-body factors are not fundamental — they are simply transmission coefficients giving the probability that modes which escape from the horizon make it out to null infinity. Past horizons are specific to eternal black holes and simply not relevant for astrophysical black holes. Similarly because the entropy-area law is tied to the validity of the Einstein equations we shall seek to avoid any appeal to this property. Since much of what I have to say is well-known to experts in the field (at least within the general relativity community) I will place a premium on clarity and simplicity.

## 2 Metric: Painlevé–Gullstrand form

In general relativity any spherically symmetric geometry, static or not, can by suitable choice of coordinates locally be put in the form

$$\mathrm{d}s^2=c(r,t)^2\mathrm{d}t^2+(\mathrm{d}rv(r,t)\mathrm{d}t)^2+r^2[\mathrm{d}\theta ^2+\mathrm{sin}^2\theta \mathrm{d}\varphi ^2].$$
(1)
Equivalently

$$\mathrm{d}s^2=[c(r,t)^2v(r,t)^2]\mathrm{d}t^22v(r,t)\mathrm{d}r\mathrm{d}t+\mathrm{d}r^2+r^2[\mathrm{d}\theta ^2+\mathrm{sin}^2\theta \mathrm{d}\varphi ^2].$$
(2)
These are so-called Painlevé–Gullstrand coordinates, a relatively obscure coordinate system currently enjoying a resurgence. (See  for a geometric discussion.) A nice feature of these coordinates is that the metric is nonsingular at the horizon. Though the basic physics is coordinate independent, we shall see that these coordinates simplify computations considerably.

In the acoustic geometries associated with sound propagation in a flowing fluid it is often convenient to not set up coordinates this precise manner but rather to use the natural coordinates inherited from the background Minkowski spacetime. Doing so results in a minor modification of the Painlevé–Gullstrand metric

$$\mathrm{d}s^2=\left[\frac{\rho (r,t)}{c(r,t)}\right]\left\{c(r,t)^2\mathrm{d}t^2+(\mathrm{d}rv(r,t)\mathrm{d}t)^2+r^2[\mathrm{d}\theta ^2+\mathrm{sin}^2\theta \mathrm{d}\varphi ^2]\right\},$$
(3)
where $`c(r,t)`$ is now the speed of sound, $`v(r,t)`$ is the radial velocity of the fluid, and $`\rho (r,t)`$ is its density. Phonons are then massless minimally coupled scalar fields in this acoustic geometry . Fortunately the conformal factor does not affect the surface gravity and does not affect the Hawking effect ; for simplicity I shall simply suppress the conformal factor, it can easily be reinstated if desired.

In matrix form the Painlevé–Gullstrand metric takes the quasi–ADM form

$$g_{\mu \nu }(t,\stackrel{}{x})\left[\begin{array}{ccc}(c^2v^2)& \mathrm{}& v\widehat{r}_j\\ \mathrm{}\mathrm{}\mathrm{}\mathrm{}& & \mathrm{}\mathrm{}\\ v\widehat{r}_i& \mathrm{}& \delta _{ij}\end{array}\right].$$
(4)
The inverse metric is

$$g^{\mu \nu }(t,\stackrel{}{x})\frac{1}{c^2}\left[\begin{array}{ccc}1& \mathrm{}& v\widehat{r}^j\\ \mathrm{}\mathrm{}& & \mathrm{}\mathrm{}\mathrm{}\mathrm{}\\ v\widehat{r}^i& \mathrm{}& (c^2\delta ^{ij}v^2\widehat{r}^i\widehat{r}^j)\end{array}\right].$$
(5)
Because of spherical symmetry finding the apparent horizon is particularly easy: it is located at $`c(r,t)=|v(r,t)|`$. It is easy to see that the metric is nonsingular at the apparent horizon, with $`det(g_{\mu \nu })=c(r,t)^2`$. To get a future apparent horizon, corresponding to an astrophysical black hole, we need $`v<0`$. I make no claims as to the location or even existence of any event horizon.

\[Since for static black holes the event horizon is coincident with the apparent horizon, the general feeling is that for quasi-stationary black holes the event horizon is likely to be “near” the apparent horizon. If accretion dominates over evaporation the event horizon is likely to be just outside the apparent horizon, while if evaporation dominates over accretion it is likely to be just inside the apparent horizon. There is an unpopular minority opinion that takes the view that because of Hawking evaporation a black hole will never form a true event horizon (absolute horizon). For the derivation presented herein this entire issue is simply not relevant.\]

Define a quantity:

$$g_H(t)=\frac{1}{2}\frac{\mathrm{d}[c(r,t)^2v(r,t)^2]}{\mathrm{d}r}|_H=c_H\frac{\mathrm{d}[c(r,t)|v(r,t)|]}{\mathrm{d}r}|_H.$$
(6)
If the geometry is static, this reduces to the ordinary definition of surface gravity :

$$\kappa =\frac{g_H}{c_H}.$$
(7)
If the geometry is not static this is a natural definition of the “surface gravity” of the apparent horizon.

## 3 Eikonal and WKB approximation ($`s`$ wave)

Consider a quantum field $`\varphi (r,t)`$ on this Painlevé–Gullstrand background and take the eikonal (geometric optics/acoustics) approximation for the $`s`$ wave

$$\varphi (r,t)=𝒜(r,t)\mathrm{exp}[i\phi (r,t)]=𝒜(r,t)\mathrm{exp}\left[i\left(\omega t^rk(r^{})dr^{}\right)\right]$$
(8)
whereby the field is written as a rapidly varying phase times a slowly varying envelope. The second equality above, where we have written the time dependence of the phase as $`\omega t`$, is valid provided the geometry is slowly evolving on the timescale of the wave, that is, provided $`\omega \text{max}\{|\dot{c}/c|,|\dot{v}/v|\}`$. (I take $`\omega `$ to be positive.)

In the eikonal approximation the d’Alembertian equation of motion becomes

$$g^{\mu \nu }_\mu \phi _\nu \phi +iϵ=0.$$
(9)
Note in particular that I have used the eikonal approximation to immediately impose Feynman’s “$`iϵ`$-prescription” on the field ($`ϵ`$ is real, positive, and infinitesimal). The metric signature is the general relativity standard $`+++`$, which is why the $`iϵ`$-prescription appears reversed relative to the particle physics standard. Also note that in invoking the prescription I have implicitly used the fact that the spacetime geometry is smooth, even at the horizon, so that it makes sense to both adopt an eikonal approximation and then within this framework use ordinary flat-space results. The use of the $`iϵ`$ prescription in this way can be traced back, at least, to the early paper of Damour and Ruffini . If one prefers, the $`iϵ`$-prescription can be rephrased in terms of analyticity of the fields on the complexified past light cone, $`𝒵(x)=\{z=x+iy\}`$, where $`x`$ is any point in spacetime and $`y`$ lies in the past light cone. This is the Lorentz-invariant generalization of analyticity in the lower half-plane, but for all practical purposes can be replaced by the $`iϵ`$ prescription.

From the preceding equation

$$\omega ^22v(r,t)\omega k[c(r,t)^2v(r,t)^2]k^2+iϵ=0.$$
(10)
So that

$$(\omega vk)^2=c^2k^2+iϵ.$$
(11)
Whence

$$\omega vk=\sigma (1+iϵ)ck;\sigma =\pm 1.$$
(12)
For specific real frequency $`\omega `$ this gives us the wave-vector $`k(r,t)`$ as

$$k=\frac{\omega }{\sigma (1+iϵ)c+v}=\frac{\sigma \omega }{(1+iϵ)c+\sigma v}=\frac{\sigma (1+iϵ)cv}{(1+iϵ)^2c^2v^2}\omega .$$
(13)
Note:

$$\sigma =+1\text{outgoing mode}.$$
(14)
$$\sigma =1\text{ingoing mode}.$$
(15)
Keeping track of the $`iϵ`$ is important only near the apparent horizon, where it is critical; everywhere else it can safely be set to zero.

While we do not really need to use the WKB approximation (physical acoustics/optics) it can be invoked at very little additional cost. The (approximate) conserved current is

$$J_\mu =|𝒜(r,t)|^2(\omega ,k,0,0).$$
(16)
Then

$$_\mu J^\mu =0|𝒜(r,t)|\frac{1}{r}.$$
(17)
So we can write

$$\varphi (r,t)\frac{𝒩}{\sqrt{2\omega }r}\mathrm{exp}\left[i\left(\omega t^rk(r^{})dr^{}\right)\right],$$
(18)
where $`𝒩`$ is some conveniently chosen normalization, to be discussed more fully below.

## 4 Near-horizon modes: ingoing

The ingoing wavenumber ($`\sigma =1`$) is

$$k_{\mathrm{in}}=\frac{\omega }{(1+iϵ)cv}.$$
(19)
Thus in the vicinity of the future apparent $`rr_H`$ (with $`vc`$) the ingoing wavevector is approximately

$$k_{\mathrm{in}}\frac{\omega }{2c_H}.$$
(20)
Thus the ingoing modes are approximately

$$\varphi (r,t)_{\mathrm{in}}\frac{𝒩_{\mathrm{in}}}{\sqrt{2\omega }r_H}\mathrm{exp}\left[i\omega \left\{t+\frac{rr_H}{2c_H}\right\}\right].$$
(21)
This means the phase velocity of the ingoing mode as it crosses the horizon (in coordinate distance per coordinate time) is $`2c_H`$. (Phase velocity equals group velocity because there is no dispersion.) We see that the ingoing modes contain no real surprises and are relatively uninteresting.

## 5 Near-horizon modes: outgoing

Now consider the outgoing mode $`\sigma =+1`$

$$k_{\mathrm{out}}=\frac{\omega }{(1+iϵ)c+v}.$$
(22)
In the vicinity of the future apparent $`rr_H`$ (with $`vc`$) the outgoing wavevector is approximately

$$k_{\mathrm{out}}\frac{\omega }{[g_H/c_H](rr_H)+iϵc_H}.$$
(23)
But because $`ϵ`$ is infinitesimal, and both $`g_H`$ and $`c_H`$ are by hypothesis positive, we can for all practical purposes rewrite this in terms of the “principal part” and a delta function contribution. That is, near the apparent horizon

$$k_{\mathrm{out}}\frac{c_H\omega }{g_H}\left\{\mathrm{}\left(\frac{1}{rr_H}\right)i\pi \delta (rr_H)\right\}.$$
(24)
Thus we can ignore the $`iϵ`$ unless we are actually intending to cross the apparent horizon. In particular, just outside the apparent horizon

$$^rk=^r\frac{dr^{}\omega }{c(r^{})|v(r^{})|}^r\frac{dr^{}c_H\omega }{g_H(r^{}r_H)}=\frac{c_H\omega }{g_H}\mathrm{ln}[rr_H].$$
(25)
Therefore (for $`r>r_H`$)

$`\varphi (r,t)_{\mathrm{out}}`$ $``$ $`𝒩_{\mathrm{out}}{\displaystyle \frac{\mathrm{exp}\left(\pm i\left[\frac{\omega c_H}{g_H}\right]\mathrm{ln}[rr_H]\right)}{\sqrt{2\omega }r_H}}\mathrm{exp}\left\{i\omega t\right\}`$
$``$ $`𝒩_{\mathrm{out}}{\displaystyle \frac{[rr_H]^{\pm i\omega c_H/g_H}}{\sqrt{2\omega }r_H}}\mathrm{exp}\left\{i\omega t\right\}.`$
The phase velocity of the outgoing mode as it crosses the horizon (in coordinate distance per coordinate time) is zero.

The fact that these outgoing modes have the surface gravity, $`\kappa =g_H/c_H`$, show up in such a fundamental and characteristic way is already strongly suggestive; and this is really all there is to Hawking radiation. The logarithmic phase pile-up at the horizon is characteristic of many derivations of Hawking radiation, in particular , and for many readers this will be sufficient to convince them that Hawking radiation is present under the current circumstances (slowly evolving apparent horizon without prejudice as to where the metric comes from). In fact, this calculation is the easiest and fastest way I know of to deduce the existence of the phase pileup using completely elementary methods.

To properly describe the modes that escape to infinity we should normalize on the half-line $`r>r_H`$ using the standard Klein-Gordon norm. This results in replacing $`𝒩_{\mathrm{out}}`$ with some specific normalization constant $`𝒩_{\mathrm{escape}}`$ whose precise value we do not need to know.

## 6 Straddling the horizon: the Boltzmann factor

In addition to looking at outgoing modes that escape to infinity, it is also useful to consider what happens to “outgoing” modes that straddle the apparent horizon. Inside the apparent horizon the real part of $`k_{\mathrm{out}}`$ goes negative, indicating that while the outgoing mode is trying to escape “upstream” it is actually being overcome by the flow/shift vector and swept back “downstream”. Additionally the phase picks up an imaginary contribution, due ultimately to the $`iϵ`$ prescription,

$`{\displaystyle _r_{}^{r^+}}k_{\mathrm{out}}`$ $``$ $`{\displaystyle _r_{}^{r^+}}dr^{}{\displaystyle \frac{c_H\omega }{g_H}}\left\{\mathrm{}\left({\displaystyle \frac{1}{r^{}r_H}}\right)i\pi \delta (r^{}r_H)\right\}`$ (26)
$`=`$ $`{\displaystyle \frac{c_H\omega }{g_H}}\left\{\mathrm{ln}{\displaystyle \frac{|r_+r_H|}{|r_{}r_H|}}i\pi \right\}.`$
So just inside the apparent horizon

$$\varphi (r,t)_{\mathrm{straddle}(r<r_H)}𝒩_{\mathrm{straddle}}\frac{|rr_H|^{\pm i\omega c_H/g_H}}{\sqrt{2\omega }r_H}\mathrm{exp}\left\{+\frac{\pi \omega c_H}{g_H}\right\}\mathrm{exp}\left[i\omega t\right].$$
(27)
Compare with the situation just outside the apparent horizon

$$\varphi (r,t)_{\mathrm{straddle}(r>r_H)}𝒩_{\mathrm{straddle}}\frac{|rr_H|^{\pm i\omega c_H/g_H}}{\sqrt{2\omega }r_H}\mathrm{exp}\left[i\omega t\right].$$
(28)
Note that this straddling mode has associated with it a normalization factor $`𝒩_{\mathrm{straddle}}`$, which is distinct from that of the escaping mode $`𝒩_{\mathrm{escape}}`$; the straddling mode is to be normalized on the entire half-line $`r>0`$. In terms of the Heaviside step function

$`\varphi (r,t)_{\mathrm{straddle}}`$ $``$ $`𝒩_{\mathrm{straddle}}\left[\mathrm{\Theta }(r_Hr)\mathrm{exp}\left\{+{\displaystyle \frac{\pi \omega c_H}{g_H}}\right\}+\mathrm{\Theta }(rr_H)\right]`$ (29)
$`\times {\displaystyle \frac{|rr_H|^{\pm i\omega c_H/g_H}}{\sqrt{2\omega }r_H}}\mathrm{exp}\left[i\omega t\right].`$
See Damour and Ruffini , equation (5b), for an early occurrence of a very similar statement; see also Massar and Parentani , equation (13). Note note in particular the presence of a Boltzmann-like factor

$$\mathrm{exp}\left\{+\frac{\pi \omega c_H}{g_H}\right\}$$
(30)
which (we shall soon see) corresponds to the Hawking temperature

$$kT_H=\frac{\mathrm{}g_H}{2\pi c_H}.$$
(31)
For many physicists the presence of this Boltzmann-like factor will be enough: The occurrence of Boltzmann factors of this type was the key to the Hartle–Hawking derivation , though they were working with the full propagator and dealing with past and future horizons of a maximally extended Kruskal–Szekeres eternal black hole. Damour and Ruffini  demonstrated the existence of similar Boltzmann factors for mode functions evaluated at the future Killing horizon of a Kerr–Newman black hole (dispensing with the past horizon entirely). In the present situation, the same Boltzmann factor is seen to arise for any slowly evolving apparent horizon.

The imaginary contribution to the integrated wavenumber is, in slightly disguised form, equivalent to the imaginary contribution to the action that occurs in the Parikh–Wilczek approach ,<sup>1</sup><sup>1</sup>1Note that Parikh and Wilczek have subsequently and implicitly made use of the Einstein equations at the stage when they then relate the emission process to the entropy change. This comment also applies to the Massar–Parentani approach. and is also related but not identical to the imaginary contribution arising from complex paths in the approach of Padmanabhan et al .

But it is possible to quite easily do a lot more: While I have so far carefully not specified specific values for $`𝒩_{\mathrm{escape}}`$ and $`𝒩_{\mathrm{straddle}}`$, the relationship between them is very simple. Since the straddling mode is to be normalized on $`(0,+\mathrm{})`$, \[which we actually approximate by the full line $`(\mathrm{},+\mathrm{})`$\], while the escaping mode is only normalized on the half line $`(r_H,\mathrm{})`$, we have

$$|𝒩_{\mathrm{straddle}}|^2\left[\mathrm{exp}\left\{+\frac{2\pi \omega c_H}{g_H}\right\}1\right]=|𝒩_{\mathrm{escape}}|^2$$

So the relative normalization is

$$\left|\frac{𝒩_{\mathrm{straddle}}}{𝒩_{\mathrm{escape}}}\right|^2=\frac{1}{\mathrm{exp}\left\{+\frac{2\pi ;\omega c_H}{g_H}\right\}1}.$$

It is this Planckian form of the normalization ratio that then leads to a Planckian distribution for the outgoing flux. (The straddling mode contains a Planckian distribution of escaping modes.)

Note the physics assumption hidden here: one is assuming that the quantum vacuum state is that corresponding to $`\varphi _{\mathrm{straddle}}`$. That is, freely falling observers (who get to see both sides of the horizon) should not see any peculiarities as one crosses the horizon. Picking the quantum vacuum corresponding to $`\varphi _{\mathrm{straddle}}`$ implies choosing the Unruh vacuum — and when we look at this vacuum state far from the horizon we see the Planckian flux of outgoing particles.

(If for whatever reason you don’t like this normalization calculation or the related thermodynamic arguments you can alternatively use the phase pileup property directly to perform a Bogolubov coefficient calculation in the style of ; all roads lead to Rome.)

## 7 Beyond $`s`$ wave

What happens if we go beyond the $`s`$ wave? There is now some momentum transverse to the apparent horizon so that

$$_\mu \phi =(\omega ,k,k_{}).$$
(32)
If we resolve the field in terms of partial waves

$$k_{}^2=\frac{\mathrm{}(\mathrm{}+1)}{r^2}.$$
(33)
Then in the eikonal approximation (I now suppress the $`iϵ`$, it has done its job and would now only serve to clutter the formulae)

$$\omega ^2+2v(r,t)\omega k+[c(r,t)^2v(r,t)^2]k^2+c(r,t)^2k_{}^2=0.$$
(34)
That is

$$(\omega vk)^2=c^2k^2+c^2k_{}^2.$$
(35)
This is a quadratic for $`k`$ as a function of $`\omega `$ and $`k_{}`$:

$$k=\frac{\sigma \sqrt{c^2\omega ^2(c^2v^2)c^2k_{}^2}v\omega }{c^2v^2}.$$
(36)
For ingoing modes near the apparent horizon one must evaluate using L’Hôpital’s rule:

$$k_{\mathrm{in}}\frac{\omega ^2c^2k_{}^2}{2c_H\omega }.$$
(37)
So the ingoing modes *do* depend on $`k_{}`$.

$$\varphi (r,t)_{\mathrm{in}}\frac{𝒩_{\mathrm{in}}}{\sqrt{2\omega }r_H}\mathrm{exp}\left[i\omega \left\{t+\frac{(rr_H)[\omega ^2c^2k_{}^2]}{2c_H\omega }\right\}\right].$$
(38)
But that does not matter: The ingoing modes are not the relevant ones. For the outgoing modes, near the apparent horizon we see

$$k_{\mathrm{out}}\frac{c_H\omega }{g_H(rr_H)}.$$
(39)
That is, for the outgoing modes:
—The near-horizon asymptotic behaviour is *independent* of $`k_{}`$.
—The phase pile-up is *independent* of $`k_{}`$.
—Continuation of the outgoing modes across horizon is *independent* of $`k_{}`$.
—The Hawking temperature *independent* of $`k_{}`$.
—This behaviour is universal for all partial waves.
—Adding a mass term corresponds to:

$$c^2k_{}^2c^2k_{}^2+(mc^2/\mathrm{})^2.$$
(40)
In view of the above we see that the behaviour of the outgoing modes near the horizon is also completely independent of the mass and transverse momentum. Consequently the *same* universal Hawking temperature applies to all masses and all partial waves. (Again, this is the easiest way I know of to convince oneself by elementary means that restricting attention to the $`s`$-wave captures almost all the essential physics of Hawking radiation.)

Note that while the Hawking temperature is completely independent of both angular momentum and mass, the grey-body factors are another matter: they do depend on both angular momentum and mass and are responsible for effectively cutting off the higher angular momentum modes.

## 8 Essential features

The only real physics input has been basic quantum physics plus the existence of a Lorentzian metric with:
— an apparent horizon;
— non-zero $`g_H`$;
— slow evolution.
The need for slow evolution of the geometry is hidden back in the approximation used to write the modes as $`\mathrm{exp}(\pm i\omega t)`$ times a position-dependent factor. This makes sense only if the geometry is quasi-static on the timescale set by $`\omega `$. So for consistency we should only trust the Boltzmann factor, and the Planckian nature of the Hawking radiation, for frequencies greater than $`\text{max}\{|\dot{c}/c|,|\dot{v}/v|\}`$. In particular, in order for the peak in the Planck spectrum to be meaningful we require

$$\frac{kT_H}{\mathrm{}}\omega _{\mathrm{peak}}\text{max}\{|\dot{c}/c|,|\dot{v}/v|\}.$$
(41)
In particular we require

$$\frac{\mathrm{d}[c(r,t)|v(r,t)|]}{\mathrm{d}r}|_H\frac{\dot{c}_H}{c_H}.$$
(42)
Near the horizon spatial gradients should dominate over temporal gradients. In particular the closer the black hole is to extremality the slower it is permitted to evolve if there is to be any hope for even a small segment of quasi-thermal spectrum.

That’s it. It is truly remarkable how basic and primitive the Hawking radiation phenomenon is, and how few physical assumptions are really necessary.

## 9 Discussion

Can the essential conditions for the Hawking effect be further relaxed?

One obvious question is the use of spherical symmetry, which precludes direct application of the current approach to Kerr and other rotating black holes. This is a technical problem, not a fundamental problem, and working in axisymmetric geometries will be do-able but somewhat messier. (For Kerr–Newman black holes the Damour–Ruffini analysis can be adapted to this end .) A tricky point for general analog model geometries is that without the Einstein equations, and something like the dominant energy condition, there is no longer any reason to believe in the zero’th law: the surface gravity and Hawking temperature can then in principle vary over different parts of the apparent horizon; these complications were suppressed in the current article via the simple expedient of enforcing spherical symmetry. (For Killing horizons there are derivations of the zero’th law that do not depend on the Einstein equations , but such considerations lose their force once the horizon becomes time dependent.)

Secondly, there are simple linguistic issues of definition: How far can we push the Hawking effect before we should give it another name? As argued in this article, based on the physics there is a very good case for keeping the name the same for the effect in arbitrary “effective geometries”, no matter how derived. Some would even argue that the Hawking effect and Unruh effects are fundamentally identical; I prefer to view then as distinct, possibly as two sides of the same coin — the response of the quantum vacuum to externally imposed conditions.

Thirdly: What is the energy source for the Hawking radiation? The infalling particles have negative energy as seen from outside the apparent horizon. In the case of general relativity black holes the infalling particles serve to reduce the total mass-energy of the central object, and it is ultimately the total mass-energy of the black hole that provides the energy emitted in the hawking flux. In the case of an acoustic black hole the infalling negative energy phonon steals kinetic energy from the fluid flow used to generate the acoustic geometry. For effective geometries associated with “slow light” the electromagnetic control field, used to generate EIT (electromagnetically induced transparency) and so reduce the group velocity, provides energy to the system which is available to ultimately be converted to Hawking-like photons. The general message is that the Hawking effect steals energy from whatever process is used to set up the effective geometry in question.

Finally, since this point still seems to cause much confusion, I should make the explicit comment that

* Hawking radiation is not a test of quantum gravity.

Instead, searching for Hawking radiation is a test of the general principles of quantum field theory in curved spacetimes. As such it is an ingredient useful for testing semiclassical quantum gravity, though it does not necessarily probe quantum gravity itself. In particular, all the proposed experimental tests of Hawking radiation via “analog models” will only probe kinematic aspects of black hole physics. To start to address the dynamics of general relativity black holes one needs the Einstein equations (or some approximation thereto), in which case one can begin to discuss Bekenstein entropy (or some approximation thereto). This requires a whole extra layer of theoretical superstructure, and a key point of this paper is that it is important and useful to keep these notions logically distinct.

###### Acknowledgments.

I acknowledge comments and feedback from the participants at the Black Holes III conference in Kananaskis, Canada. I particularly wish to thank Bei-Lok Hu, Ted Jacobson, Renaud Parentani, and Amanda Peet for their comments.


text:
# Extracting the Flavour Structure of a Baryon-Number R-parity Violating Coupling at the LHC

## 1 Introduction

In supersymmetry (SUSY) models where R-parity, $`\mathrm{R}_\mathrm{p}=(1)^{3\mathrm{B}+\mathrm{L}+2\mathrm{S}}`$, is violated (RPV models) the lightest supersymmetric particle (LSP) can decay to Standard Model particles. For the baryon-number violating RPV coupling, $`\lambda ^{\prime \prime }`$, the LSP, assumed to be the $`\stackrel{~}{\chi }_1^0`$, decays via a virtual squark into three quarks (or antiquarks). In previous work  we suggested a method by which the masses of several of the SUSY particles can be measured in such scenarios for a wide range of $`\lambda ^{\prime \prime }`$.

If R-parity violating couplings are observed, it seems probable that they come from an underlying theory which also determines the fermion mass spectrum. This idea has been explored  in the context of gauged $`U(1)_X`$ models of flavour, in which a hierarchical pattern of RPV couplings is generated by the same symmetry from which the fermion mass spectrum is derived. Thus the determination of the flavour structure of an RPV coupling could be a vital clue in the resolution of the flavour problem.

In this paper we investigate the flavour content of the quark jets from $`\stackrel{~}{\chi }_1^0`$ decays. The aim is to identify which $`\lambda _{ijk}^{\prime \prime }`$ coupling is dominant in the neutralino decay. The index $`i`$ is the generation number of an up-type quark, and $`j`$ and $`k`$ correspond to down-type quarks. The coupling is antisymmetric in $`j`$ and $`k`$, so there are nine possible non-zero elements.

Three of the couplings, $`\lambda _{\mathrm{t}jk}^{\prime \prime }`$, produce very different signatures because they allow the $`\stackrel{~}{\chi }_1^0`$ to decay via the top quark. If the neutralino is lighter (heavier) than the top quark then the decay is propagator (phase-space) suppressed. For the mSUGRA point we investigate (which is described in Section 2) the branching ratio through the top coupling is a factor of $`4.2\times 10^6`$ smaller than through a non-top coupling of the same size.

For a top mode to dominate the decay, a $`\lambda ^{\prime \prime }`$ coupling to top would need to be a factor of at least 2000 greater than any of the non-top couplings. In that case the lightest neutralino would typically be long-lived (Figure 1) and for $`\lambda _{\mathrm{t}jk}^{\prime \prime }<\mathrm{\hspace{0.17em}\hspace{0.17em}10}^2`$ nearly all $`\stackrel{~}{\chi }_1^0`$ s decay outside the detector. Determination of these couplings could then only be achieved by searching for the rare decay of $`\stackrel{~}{\chi }_1^0`$s in the active volume, while sparticle mass measurements could be made in the same way as for the $`\mathrm{R}_\mathrm{p}`$ conserving case.

We investigate the other six couplings – $`\lambda _{\mathrm{uds}}^{\prime \prime }`$, $`\lambda _{\mathrm{udb}}^{\prime \prime }`$, $`\lambda _{\mathrm{usb}}^{\prime \prime }`$, $`\lambda _{\mathrm{cds}}^{\prime \prime }`$, $`\lambda _{\mathrm{cdb}}^{\prime \prime }`$ and $`\lambda _{\mathrm{csb}}^{\prime \prime }`$ – assuming that one is significantly larger than the others. Since the branching ratio is proportional to the square of the coupling, it is sufficient for this analysis that all the subdominant couplings are 4.5 times smaller than the dominant coupling. Then just 5% of all $`\stackrel{~}{\chi }_1^0`$s will decay though the subdominant channel, and 90% of events, each with two $`\stackrel{~}{\chi }_1^0`$s, will contain only dominant-coupling decays.

## 2 Event Simulation and Selection

HERWIG 6.2 is used as the event generator, and the official ATLAS simulation program ATLFAST 2.50 , is used to simulate the performance of the ATLAS detector. We use as a test case the mSUGRA point $`m_0=100`$ GeV, $`m_{1/2}=300`$ GeV, $`A_0=300`$ GeV, $`\mathrm{tan}\beta =10`$ and $`\mathrm{sgn}\mu +`$, with one of $`\lambda _{ijk}^{\prime \prime }=5\times 10^3`$. This is in the middle of the range of couplings $`10^5<\lambda ^{\prime \prime }<\mathrm{\hspace{0.17em}\hspace{0.17em}1}`$, for which we can expect double $`\stackrel{~}{\chi }_1^0`$ decay within the detector, and is approaching indirect limits set from nuclear decays for $`\lambda _{\mathrm{uds}}^{\prime \prime }`$ and $`\lambda _{\mathrm{udb}}^{\prime \prime }`$ . A selection of the sparticle masses for this point is shown in Table 1. The applicability of our analysis to SUGRA models with other parameters is discussed in our previous work . The method of the event selection, which is also detailed in that paper, is summarised below.

Signal events will contain three jets from each of the two $`\stackrel{~}{\chi }_1^0`$s. Since the SUSY cross section is dominated by squark and gluino production, there are usually at least two other jets from the $`\mathrm{R}_\mathrm{p}`$ conserving decays of squarks and gluinos. Requiring the presence of leptons, which can originate in the decay chain $`\stackrel{~}{\chi }_2^0\stackrel{~}{l}_R^{}\mathrm{}^\pm \stackrel{~}{\chi }_1^0\mathrm{}^{}\mathrm{}^\pm `$, helps to decrease the Standard Model (SM) background. Events are selected with between eight and ten jets, and with two opposite-sign, same family isolated leptons (e or $`\mu `$).

Cuts on event-shape variables are applied, reducing the SM background to less than 10% of the SUSY signal, and kinematic cuts are applied to preferentially select jets from neutralino decays. All possible three-jet combinations are inspected and their invariant masses calculated. Since we expect the event to contain two LSP decays we select two of these three-jet combinations with reconstructed masses within 20 GeV of each other as $`\stackrel{~}{\chi }_1^0`$ candidates. To reconstruct the $`\stackrel{~}{\chi }_2^0`$ we pick the $`\stackrel{~}{\chi }_1^0`$ candidate closest in $`\eta \varphi `$ to either of the leptons from the lepton pair, and find the total invariant mass of the three jets plus those two leptons.

A 2-dimensional Gaussian is fitted to the $`\stackrel{~}{\chi }_1^0`$ and $`\stackrel{~}{\chi }_2^0`$ masses, and combinations are selected within $`2\times \sigma `$ of that peak. This selection removes much of the combinatoric background, so that the jet multiplicity cut can be relaxed to $`8N_{\mathrm{jet}}12`$. We now have a rather clean sample of events in which we can identify six jets with the two neutralino decays. All $`N_{\mathrm{comb}}`$ jet combinations which pass the above cuts for any event are accepted and are given weight $`1/N_{\mathrm{comb}}`$.

## 3 Vertex Tagging

The lifetimes of hadrons containing b- and c-quarks give rise to displaced vertices which can be reconstructed from charged tracks in the inner detector. In general the lifetimes of charmed hadrons are shorter than for hadrons containing a b-quark, for example $`c\tau _{D^0}=124\mu \mathrm{m}`$ compared to $`c\tau _{B^0}=468\mu \mathrm{m}`$. This allows statistical separation of c-quark jets from b-quark jets, and of c- and b-jets from light quark (u, d, s) and gluon jets. Strange hadrons have longer lifetimes with $`c\tau `$ of the order of tens of centimeters. Since they decay to a small number of particles it is difficult to reconstruct a secondary vertex, and so strange hadrons cannot be tagged within jets in the LHC environment.

The vertex tagging performance of the ATLAS Inner Detector was simulated in . In that study a likelihood ratio method was applied to the transverse impact parameter $`(d_0)`$ of the tracks within the jet cone. The resultant b-tagging efficiency and the rejection rates for c- and other jets were parameterised as functions of the transverse momentum and pseudorapidity of the reconstructed jet and were implemented in ATLFAST. The rejection factors, averaged over all directions and transverse momenta, are summarised in Table 2. In order to select a higher purity sample, we tolerate a rather low b-tagging efficiency (33%).

The proportions of events with 0, 1 and 2 tagged jets are shown in Figure 2, for the six different $`\lambda ^{\prime \prime }`$ couplings. As we expect, the vertex tagging rate is greatest when the $`\stackrel{~}{\chi }_1^0`$ decay products include both c and b quarks, and smallest for couplings which produce light-quark daughters. As we cannot distinguish d- from s-jets, there is no discrimination between $`\lambda _{i\mathrm{db}}^{\prime \prime }`$ and $`\lambda _{i\mathrm{sb}}^{\prime \prime }`$ (for $`i`$=u,c).

In these simulations, with $`\lambda ^{\prime \prime }=5\times 10^3`$ the $`\stackrel{~}{\chi }_1^0`$ lifetime of $`10^{14}`$ seconds corresponds to a typical decay length $`c\tau `$ of $`3\mu \mathrm{m}`$. This is short in comparison with typical b- and c-jet vertex displacements so it should not severely affect vertex reconstruction. However, the lifetime is inversely proportional to the square of the coupling strength. When $`c\tau `$ becomes bigger than about 3 cm (for $`\lambda ^{\prime \prime }<\mathrm{\hspace{0.17em}\hspace{0.17em}5}\times 10^5`$) special vertex reconstruction would be required. Vertex tagging becomes almost impossible when neutralinos travel more than about 30 cm in the transverse direction since their daughters will produce hits only in the outermost silicon layers.

## 4 Flavour Discrimination from Muons

The weak decay of hadrons containing heavy quarks can produce muons, which will generally lie within the associated quark jet. The muons will pass though the calorimeter and so can be measured by the muon detector even if they lie inside the jet cone. The frequency with which these “non-isolated” muons occur in signal events can be used as an additional discriminator between quark flavours.

The number of non-isolated muons per event in the six signal jets is plotted in Figure 3. Since bottom mesons produce more muons \[$`B(B^0\mu ^+\nu _\mu \mathrm{X})=10\%`$\] than charmed mesons \[$`B(D^0\mu ^+\nu _\mu \mathrm{X})=7\%`$\] , we can use the number of muons in signal jets to statistically separate b- from c-jets.

## 5 Statistical Significance

The confidence with which we can identify the dominant coupling was explored for all pairs of couplings for both the vertex tagging and muon rates (Figures 2 and 3 respectively). The variable,

$`\chi ^2={\displaystyle \underset{i}{}}{\displaystyle \frac{(x_iy_i)^2}{\sigma _x^2+\sigma _y^2}}`$
was calculated, where $`x_i`$ and $`y_i`$ are the fractions of events with $`i`$ muons (or vertex tags) for the test pair of couplings and the $`\sigma ^2`$s are their variances assuming Poisson-distributed numbers of events. Bins where one or other distribution contained fewer than five events were excluded.

The calculated $`\chi ^2`$ values for both the muon and vertex-tagging plots are shown in Table 3. If we use only the muon information then all couplings are distinguishable at 90% except for the ambiguity between d- and s-jets. In general the vertexing rate gives better discrimination, but would be difficult for a long-lived $`\stackrel{~}{\chi }_1^0`$ for the reasons discussed in Section 3. Combining the results from both analyses gives separation at better than 3.5 $`\sigma `$ in all cases, again with the exception that it is not possible to distinguish down from strange quarks.

If the branching ratios of the neutralino through two or more couplings were comparable, then as the information in the distributions is rather degenerate, definitive identification of the couplings becomes difficult. However the method could then be used to constrain the possible values of the couplings.

## 6 Conclusions

If baryon-number violating couplings to quarks other than top are $`>\mathrm{\hspace{0.17em}\hspace{0.17em}5}\times 10^5`$ then $`\stackrel{~}{\chi }_1^0`$s will typically decay within the tracking volume of an LHC detector. For these couplings, our simulations show that displaced secondary vertices and muons from heavy quark decays allow statistical separation of b- from c-quark jets, and b- and c-jets from light-quark jets. A dominant RPV coupling can then be identified at better than 3.5 $`\sigma `$, except for an ambiguity caused by the inability to distinguish strange from down quarks.

We have demonstrated that if the MSSM breaks baryon number, not only can the LHC detect and measure the masses of sparticles, but it can even distinguish the flavour structure of the RPV coupling.

###### Acknowledgments.

We would like to thank F. Paige for suggesting this extension of our previous work. We thank H. Dreiner, members of the Cambridge SUSY working group, and the ATLAS collaboration for helpful discussions. We have made use of the physics analysis framework and tools which are the result of collaboration-wide efforts. This work was partly funded by PPARC.


text:
# 1 Parton luminosities as a function of 𝜏=𝑠̂/𝑠 at the LHC. The inset plot shows the luminosities at the Tevatron as a function of the partonic centre-of-mass energy √𝑠̂.

Although the Higgs boson is a crucial ingredient of the Standard Model (SM), it has so far eluded experimental discovery. LEP results  imply a lower limit of $`M_H>113.5`$ GeV (at $`95\%`$ CL) on the mass $`M_H`$ of the SM Higgs boson, favour a light Higgs boson ($`M_H<200`$ GeV) and suggest the observation of signal events at $`M_H115`$ GeV.

After the end of the LEP program, the Higgs search will be carried out at hadron colliders, where the dominant production mechanism is gluon–gluon fusion through a heavy-quark loop. At the LHC , $`gg`$ fusion overwhelms the other production channels in the case of a light Higgs boson, and it still provides $`50\%`$ of the total production rate at $`M_H1`$ TeV. At the Tevatron , $`gg`$ fusion gives $`65\%`$ of the total cross section in the mass range $`M_H=100`$–200 GeV, although this production mechanism leads to an important discovery mode only at $`M_H>135`$ GeV, when the decay channel $`HWW^{}`$ opens up (if $`M_H<135`$ GeV, the decay channel $`Hb\overline{b}`$ suffers from a huge QCD background, and the decay rate $`H\gamma \gamma `$ is too low to be observed).

The NLO QCD corrections to Higgs boson production through $`gg`$ fusion give a large effect . Since approximate evaluations  of higher-order terms suggest that their effect can still be sizeable, the computation of the NNLO corrections is certainly important.

A first step in this direction has been performed by two groups , who have computed the soft and virtual contributions to the NNLO partonic cross section $`\widehat{\sigma }(ggH+X)`$ in the large-$`m_{\mathrm{top}}`$ approximation. Our calculation  was done by combining the recent result  for the two-loop amplitude $`ggH`$ with the soft factorization formulae for tree-level  ($`ggHgg,Hq\overline{q}`$) and one-loop  ($`ggHg`$) amplitudes. The independent calculation of Ref.  uses a different method, and the results for the partonic cross section fully agree.

In Ref. , we have combined the NNLO partonic cross section with the recent MRST2000 set of parton distributions , which includes (approximate) NNLO densities, thus providing a first consistent estimate of the NNLO QCD corrections to the hadronic cross section. In the following, we briefly summarize our main results.

The soft and virtual terms give the most important contributions to the partonic cross section when it is computed at a value of the partonic centre-of-mass energy $`\sqrt{\widehat{s}}`$ that is close to the Higgs boson mass $`M_H`$. The hadronic centre-of-mass energy $`\sqrt{s}`$ is related to the partonic one as $`\widehat{s}=\tau s=x_1x_2s`$, where $`x_1`$ and $`x_2`$ are the momentum fractions of the two partons that initiate the partonic subprocess. The density of these partons is controlled by the parton distributions of the colliding hadrons. Since the parton distributions are strongly suppressed at large values of the momentum fractions $`x_1,x_2`$, $`\widehat{s}`$ is typically much smaller than $`s`$ and, therefore, the soft and virtual terms can approximate the complete result also when the hadronic energy $`\sqrt{s}`$ is not very close to $`M_H`$. At fixed $`M_H`$, the quantitative reliability of the approximation depends on the value of $`\sqrt{s}`$ and on the actual value of the parton luminosity. The $`gg`$ and $`(q+\overline{q})g`$ luminosities are shown in Fig. 1. At the LHC, where we use the reference value $`M_H=115`$ GeV, the parton luminosity decreases by about two orders of magnitude when $`\sqrt{\widehat{s}}`$ increases from $`M_H`$ to $`3M_H`$ (i.e. when $`\tau `$ increases from $`10^4`$ to $`10^3`$). At the Tevatron (see the inset plot), the parton luminosity decreases by almost two orders of magnitude when $`\sqrt{\widehat{s}}`$ increases from $`\sqrt{\widehat{s}}=M_H=150`$ GeV to $`300`$ GeV. This observation suggests that the soft and virtual terms can give a good approximation at the LHC and a (slightly) better approximation at the Tevatron Run II.

In Ref. , we introduced two approximations, named soft–virtual (SV) and soft–virtual–collinear (SVC), of the partonic cross section. In the SV approximation only the contributions of soft and virtual origin are taken into account. The SVC approximation extends the SV approximation by including the next-to-dominant contribution (which has a collinear origin ) in the expansion around the region where $`\sqrt{\widehat{s}}=M_H`$.

In the following we present numerical results for the hadronic cross section evaluated up to LO, NLO and NNLO. At the NLO we compare the exact result with those of the SV and SVC approximations. At the NNLO, we add the exact NLO result to our SV and SVC approximations of the NNLO corrections. All the results are presented in term of K-factors, which are defined as the ratio of the hadronic cross section over its LO value. The latter is computed by fixing the factorization and renormalization scales $`\mu _F`$ and $`\mu _R`$ at the central value $`\mu _F=\mu _R=M_H`$. As a reference, we use the parton distributions of the MRST2000 set, with parton densities and $`\alpha _\mathrm{S}`$ evaluated at each corresponding order.

At the LHC, we find that the results do not strongly depend on the choice of parton distributions: very similar results are obtained by using the CTEQ5 set , whereas larger ($`10\%`$) deviations appear by using the GRV98 set . At the Tevatron Run II, the dependence on the parton distributions is stronger than at the LHC. The CTEQ5 set gives a LO cross section at $`M_H=100`$ GeV ($`M_H=200`$ GeV) that is $`10\%`$ ($`30\%`$) smaller than the value obtained by using the MRST2000 set. Of course, this dependence affects the K-factors through their normalization with respect to the LO cross section.

The K-factors at the LHC and at the Tevatron are shown in Fig. 2. The bands are obtained by independently varying $`\mu _F`$ and $`\mu _R`$ in the range $`M_H/2\mu _F,\mu _R2M_H`$. The results show a reduction of the scale dependence at NNLO (from about $`\pm 20\%`$ at full NLO to about $`\pm 10\%`$ and $`\pm 15\%`$ at NNLO-SV and NNLO-SVC, respectively).

From the plots on the left-hand side of Fig. 2, we see that at NLO the SV approximation tends to underestimate the exact result, whereas the SVC approximation slightly overestimates it. Therefore, it is reasonable to expect the exact NNLO K-factor to lie inside the corresponding SV and SVC bands. In the case of the production of a light Higgs boson ($`M_H=100`$–200 GeV) at the LHC (Tevatron), this expectation corresponds to $`K2.2`$$`2.4`$ ($`K3`$), i.e. to an enhancement of the NLO cross section by about 10–25$`\%`$ ($`50\%`$).

Note that, at fixed $`M_H`$, Higgs boson production at the Tevatron is closer to threshold than at the LHC. This implies (see Fig. 2) that the approximation of the fixed-order results in terms of the corresponding SV and SVC contributions works better at the Tevatron than at the LHC. This also implies that, order by order in perturbation theory, QCD radiative corrections are larger at the Tevatron than at the LHC. In particular, QCD contributions beyond NNLO, due to multiple soft-gluon emission, can still be numerically relevant to Higgs boson production at the Tevatron.

Acknowledgments.
This work was supported in part by the EU Fourth Framework Programme “Training and Mobility of Researchers”, Network “Quantum Chromodynamics and the Deep Structure of Elementary Particles”, contract FMRX–CT98–0194 (DG 12 – MIHT).