text:
# LiV2O4: frustration induced heavy fermion metal

## Abstract

We propose a two-stage spin-quenching scenario for the unusual heavy fermion state realized in the mixed valent metal LiV2O4. In this theory, local valence fluctuations are responsible for the formation of partially quenched, spin $`\frac{1}{2}`$ moments below room temperaturekondo . Frustration of the intersite spin couplings then drives the system to produce the heavy Fermi liquid seen at low temperatures.takagi  The anomalous resistivity and the sign change of the Hall constant can be understood naturally within this model, which also predicts a unique symmetry for the heavy quasiparticle bands that that may be observed in de Haas-van Alphen experiments.

Magnetic frustration has attracted increasing attention in recent years as a new tool to suppress antiferromagnetism and generate new types of electronic behavior. One exciting possibility is the development of frustration induced heavy electron behavior in d-electron systems where the Kondo temperature is generally too small to overcome magnetic order without frustration. The recognition of heavy fermion behavior in LiV2O4kondoprl  was an important first step in this direction. In LiV2O4, the frustrated lattice remains undistorted to the lowest temperatures measured,xray  The heavy electron state which then develops is unusual in many respects: it displays a monotonically increasing resistivitytakagi , a pressure-driven metal-insulator transitiontakag  and field independent heat capacity up to 30Ttakag .

LiV2O4 forms a spinel structure. The magnetic vanadium atoms are homogeneously mixed valent forming a lattice of undistorted corner-shared tetrahedra. The spin at each site fluctuates between S = 1 (3d2) and S = $`\frac{1}{2}`$ (3$`d^1`$) with a formal valence of $`3.5`$. Antiferromagnetic couplings between these spinsmuhtar ; lee  imply that the lattice is a highly frustrated metal. There is a strong on-site Hund’s coupling between spins of order 1 eVanisimov , which is manifested in the photoemissionfujimori  and the high temperature magnetic susceptibilitymuhtar . At room temperature the magnetic suceptibility and specific heat indicate the presence of spin $`\frac{1}{2}`$ local moments at every site.kondo ; takagi . Below 4K, the moments quench into a heavy Fermi liquid.kondo

There have been several attempts to explain this physics. One class of model sees the frustration as secondary, the essential physics being due to mixed valencevarma , exhaustionanisimov , 1D chainsfujimoto  or a two-band Hubbard modelkusunose . Another argues that frustration is essential, leading to a tetrahedron rulefulde , a frustrated Kondo latticeburdin , antiferromagnetically coupled spins Hund’s coupled to conduction electronslacroix  or quantum criticalityshannon . Here we present a simple two-band model which we believe captures the essential physics: a two stage spin quenching that could not occur in the absence of frustration; and show how our results can be extended to the pyrochlore structure.

The high temperature physics of LiV2O4 can be understood as an S=1 Mott insulator doped with holes on a pyrochlore lattice. The strong Hund’s coupling leads us to consider only spin-1 and spin-$`\frac{1}{2}`$ sites coupled antiferromagnetically over the lattice via direct exchangesingh ,

$`H`$ $`=`$ $`t{\displaystyle \underset{(i,j)}{}}(X_{0\sigma }(i)X_{\sigma 0}(j)+h.c.)+{\displaystyle \underset{(i,j)}{}}Js_is_j`$ (1)
where $`X_{0\sigma }=|\frac{1}{2},\sigma ><1,\sigma ^{}|=\chi ^{}b_\sigma `$ is a Hubbard operator allowing the charged holes to move in the spin background via a hopping matrix element $`t`$, subject to the constraint $`b_\sigma ^{}b_\sigma +\chi ^{}\chi =2`$. $`J`$ is an antiferromagentic coupling between nearest neighbor spins. The local physics of this model is that of a mixed valent impurity undergoing charge fluctuations $`d^2d^1+e^{}`$ via hybridization with the surrounding bath of holes, as might be formalized in a dynamical mean-field description of the above model. The effective impurity Hamiltonian replacing the first term in $`H`$ would then be

$$H_1=\underset{k\sigma }{}ϵ_kc_{k\sigma }^{}c_{k\sigma }+\underset{k\sigma }{}t(X_{\sigma 0}c_{k\sigma }^{}+h.c.).$$
(2)
As the frustration of the pyrochlore lattice prevents magnetic ordering from occuring, this term dominates the high temperature physics. From strong coupling argumentskrishna-murthy , we know that a Fermi sea of electrons forms below $`T_1`$, quenching the local moment character on each site to S=$`\frac{1}{2}`$. This provides a natural explanation for the behavior seen below room temperature. On the lattice, this fluid can be visualized as a smooth Fermi sea, from which the residual spin $`\frac{1}{2}`$ moments protrude like spines on the hide of a porcupine. (Fig. 1.) The original antiferromagnetic interaction between nearest neighbor sites now couples the electron sea to the residual spin $`\frac{1}{2}`$ moments. The properties of this “porcupine” Fermi sea are the topic of the remainder of this paper. The question is – what is the nature of the intrasite screening process that removes the free spins (spines) at low temperatures?

At low temperatures, $`J`$ becomes an effective Kondo coupling between the localized spins and the newly freed electrons in the Fermi sea. The locally symmetric spin state required to satisfy the strong Hund’s coupling means that spatial overlap of the realized Wannier states is forbidden as the heavy quasiparticle forms at low temperatures. Thus, even after the complete break-down of any orbital picture, as the mixed valent spin quenching of the upper level occurs in the totally symmetric channel the second stage spin quenching can occur only in an orthogonal channel. This is the analogue of the atomic Hund’s selection rule for itinerant states. Thus the conditions

$`\{\begin{array}{cc}\text{no overlap}\hfill & \\ \text{no phase transition (equivalent sites)}\hfill & \\ \text{no time reversal symmetry breaking}\hfill & \end{array}`$
arise, uniquely specifying the second channel (Fig.2). The unusual symmetry of the composite quasi-particles resulting from this second spin-quenching should be observable by de Haas - van Alphen experiments.

We now treat the second stage at the mean field level keeping in mind that for LiV2O4, we really want to have a half-filled localized d-band interacting with a quarter-filled delocalized d-band. Estimates of $`J\theta _{CW}^{highT}400K`$muhtar , T$`{}_{c2}{}^{}450K`$ and $`t2eV`$singh  (band theory) have been made, and we expect that $`T_{c1}t`$ and $`T_{c2}te^{\frac{\mathrm{\#}tb^2}{J}}`$ as these are the two energy scales in the problem. It is convenient to recast the Hubbard operators into a slave boson language,

$`H`$ $`=`$ $`{\displaystyle \underset{<ij>}{}}(tX_{0\sigma }^iX_{\sigma 0}^j+J𝐬_i𝐬_j)`$ (3)
$``$ $`J_H{\displaystyle \underset{j}{}}(𝐬_j^22+{\displaystyle \frac{5}{4}}X_{00}^j)`$
where $`X_{0\sigma }=b^{}c_\sigma `$ is a Hubbard operator allowing the holes to move in the upper e$`{}_{}{}^{}{}_{g}{}^{}`$ symmetry channel, $`𝐬=\left(\begin{array}{cc}c_\alpha ^{}& d_\alpha ^{}\end{array}\right)\frac{\sigma _{\alpha \beta }}{2}\left(\begin{array}{c}c_\beta \\ d_\beta \end{array}\right)`$ is the total spin, and we need to enforce the constraints $`n_c+n_b=1`$, $`n_d=1`$ (a localized spin-$`\frac{1}{2}`$ level of $`a_{1g}`$ symmetry). The Hund’s coupling $`J_H`$ is to be taken to infinity to enforce the constraint $`S^2=2`$ at doubly occupied sites.

We now demonstrate how a mean field treatment of Eq.3 leads us to the effective Hamiltonian:

$`H_{mf}`$ $`=`$ $`{\displaystyle \underset{k\sigma }{}}\left(\begin{array}{cc}𝒞_{k\sigma }^{}& 𝒟_{k\sigma }^{}\end{array}\right)𝐡_k\left(\begin{array}{c}𝒞_{k\sigma }\\ 𝒟_{k\sigma }\end{array}\right)+2{\displaystyle \underset{\mathrm{\Gamma }}{}}{\displaystyle \frac{\overline{\mathrm{\Delta }}_\mathrm{\Gamma }\mathrm{\Delta }_\mathrm{\Gamma }}{J^\mathrm{\Gamma }}}`$
$`+`$ $`\lambda _1(|b|^21)\lambda _2`$
where

$`𝐡_k=\left(\begin{array}{cc}2t|b|^2\mathrm{\Phi }_k^{(1)}+\frac{\lambda _1\mu }{N}& \mathrm{\Delta }_\mathrm{\Gamma }\mathrm{\Phi }_k^{(\mathrm{\Gamma })}+i\xi \\ \overline{\mathrm{\Delta }}_\mathrm{\Gamma }\mathrm{\Phi }_k^{(\mathrm{\Gamma })}+i\overline{\xi }& \frac{\lambda _2}{N}\end{array}\right)`$ (5)
and $`𝒟_k^{}=\left(\begin{array}{ccc}d_{1k}& ..& d_{Nk}\end{array}\right),𝒞_k^{}=\left(\begin{array}{ccc}c_{1k}& ..& c_{Nk}\end{array}\right)`$ demark the number of bands arising from the magnetic unit cell of the lattice under consideration. $`\lambda _1`$ and $`\lambda _2`$ are Lagrange multipliers fixing the occupancies of the electron and spin liquid. $`\mathrm{\Phi }_k^{(\mathrm{\Gamma })}`$ is an N dimensional matrix representing the Wannier functions, where N is the number of atoms per unit cell and $`\mathrm{\Gamma }=1,..,M`$ specifies the number of orthogonal channels required to decouple the Heisenberg interaction between one site and its $`M`$ nearest neighbors.

The mean field kinetic energy term of $`H`$ is

$`H_{kin}={\displaystyle \underset{<ij>}{}}t(c_\sigma ^{}b)_i(b^{}c_\sigma )_j=2t|b|^2{\displaystyle \underset{k}{}}𝒞_k^{}\mathrm{\Phi }_k^{(1)}𝒞_k`$ (6)
where the bosons take an expectation value, and the fourier transform happens to give $`\mathrm{\Phi }_k^{(1)}`$, the totally symmetric Wannier state. Magnetic frustration enables us to neglect the normally dominant interaction between the residual spin $`\frac{1}{2}`$ moments, however, we must still consider the Kondo coupling between localized and delocalized bands in the expansion of the $`Js_is_j`$ term of Eq.3:

$`H_J=J{\displaystyle \underset{<ij>}{}}𝐬_i^d𝐬_j^c{\displaystyle \underset{i,\mathrm{\Gamma }}{}}{\displaystyle \frac{J^\mathrm{\Gamma }}{2}}(d_{i\alpha }^{}\psi _{i\mathrm{\Gamma }\alpha })(\psi _{i\mathrm{\Gamma }\beta }^{}d_{i\beta })`$ (7)
where the operators of the delocalized level have been expanded in terms of orthogonal linear combinations as $`\psi _{i\mathrm{\Gamma }}^{}=_j\varphi ^\mathrm{\Gamma }(ji)c_j^{}`$ and $`\varphi ^\mathrm{\Gamma }`$ is a form factor reflecting the phases at different nearest neighbor sites of i, reflecting the possible arrangements of c on the lattice. This becomes:

$`H_J={\displaystyle \underset{i,\mathrm{\Gamma }}{}}\left[(\mathrm{\Delta }_{i\mathrm{\Gamma }}\psi _{i\mathrm{\Gamma }\beta }^{}d_{i\beta }+\overline{\mathrm{\Delta }}_{i\mathrm{\Gamma }}d_{i\beta }^{}\psi _{i\mathrm{\Gamma }\beta })+2{\displaystyle \frac{\overline{\mathrm{\Delta }}_{i\mathrm{\Gamma }}\mathrm{\Delta }_{i\mathrm{\Gamma }}}{J^\mathrm{\Gamma }}}\right]`$ (8)
where we have performed a Hubbard-Stratonovich gauge transformation on the four-fermion term before fourier transforming to obtain the off-diagonal terms in Eq.(5).

To treat the term of interest from a Hund’s coupling:

$`H_H={\displaystyle \frac{J_H}{4}}(c_\alpha ^{}d_\alpha d_\beta ^{}c_\beta ){\displaystyle \frac{4}{J_H}}\overline{\xi }\xi +i\overline{\xi }d_\alpha ^{}c_\alpha +ic_\alpha ^{}d_\alpha \xi `$ (9)
a Hubbard-Stratonovich gauge field has been introduced to decouple the four-fermion termwhere . Before the decoupling, $`H_H|S=1>=0|S=1>`$ and $`H_H|S=0>=\frac{J_H}{2}|S=0>`$, so in the limit $`J_H\mathrm{}`$, the $`|S=0>`$ state will be forbidden. To write Eq.(4) we have taken the limit $`J_H\mathrm{}`$. Although this appears to be non-Hermitian, the mean field values of $`\xi `$, $`\overline{\xi }`$ (independent variables) will also be imaginary so that physical quantities are real.

The Free energy can then be expressed as

$`F`$ $`=`$ $`T{\displaystyle \underset{k,\sigma }{}}Tr[ln(1+e^{\beta 𝐡})]+2{\displaystyle \underset{\mathrm{\Gamma }}{}}{\displaystyle \frac{\overline{\mathrm{\Delta }}_\mathrm{\Gamma }\mathrm{\Delta }_\mathrm{\Gamma }}{J^\mathrm{\Gamma }}}`$ (10)
$`+`$ $`\lambda _1(|b|^21)\lambda _2`$ (11)
where $`𝐡_k`$ is the Hamiltonian expressed as a matrix as above. Since $`\overline{\xi }`$ and $`\overline{\mathrm{\Delta }}`$ both couple to $`𝒞_\alpha ^{}𝒟_\alpha `$, changing one affects the other, such that $`\frac{^2F}{\overline{\mathrm{\Delta }}\mathrm{\Delta }}\delta \overline{\mathrm{\Delta }}+\frac{^2F}{\overline{\xi }\mathrm{\Delta }}\delta \overline{\xi }=0`$. Integration over such fluctuations yields $`\frac{1}{J^{}}=\frac{1}{J}\frac{F_{\overline{\mathrm{\Delta }}\xi }F_{\overline{\xi }\mathrm{\Delta }}}{2F_{\xi \xi }}`$, which guarantees the orthogonality of the second channel. Channels with $`F_{\overline{\mathrm{\Delta }}\xi }0`$ are effectively removed from consideration. For example, on the square lattice at $`\frac{1}{4}`$-filling of the upper level, in the absence of the Hund’s term, $`\mathrm{\Phi }^1`$ would have the highest Tc2 while including the Hund’s term the calculated Tc2 for this channel turns out to be negative. That is no transition would occur if channels with F$`{}_{\overline{\mathrm{\Delta }}\xi }{}^{}=0`$ did not exist.

To solve for Tc2, it is helpful to isolate the interaction contribution to the Green’s functions for the delocalized level

$`G_c^1=(G_c^0)^1\mathrm{\Sigma }`$ (12)
where $`\mathrm{\Sigma }=\text{}=|\mathrm{\Delta }_\mathrm{\Gamma }|^2\mathrm{\Phi }_k^\mathrm{\Gamma }(\omega \lambda _2)^1\mathrm{\Phi }_k^\mathrm{\Gamma }`$ describes the interactions with the localized level. Writing the free energy contribution as $`F_{|\mathrm{\Delta }_\mathrm{\Gamma }|^2}=TTr\mathrm{ln}(G_c^1)+\frac{2|\mathrm{\Delta }_\mathrm{\Gamma }|^2}{J^\mathrm{\Gamma }}`$ we find,

$`{\displaystyle \frac{^2F}{\mathrm{\Delta }_\mathrm{\Gamma }\overline{\mathrm{\Delta }}_\mathrm{\Gamma }}}|`$ $`{}_{\mathrm{\Delta }_\mathrm{\Gamma }=i\xi =0}{}^{}={\displaystyle \frac{2}{J^\mathrm{\Gamma }}}+T{\displaystyle \underset{n,k}{}}Tr[\mathrm{\Phi }_k^\mathrm{\Gamma }G_c^0G_d^0\mathrm{\Phi }_k^\mathrm{\Gamma }]`$
$`=`$ $`{\displaystyle \underset{k}{}}{\displaystyle \underset{\alpha =1}{\overset{N}{}}}{\displaystyle \frac{f(E_{k\alpha })f(\lambda _2)}{\lambda _2E_{k\alpha }}}\eta _{k\alpha }^{}\mathrm{\Phi }_k^\mathrm{\Gamma }\mathrm{\Phi }_k^\mathrm{\Gamma }\eta _{k\alpha }+{\displaystyle \frac{2}{J^\mathrm{\Gamma }}}=0`$
where $`E_{k\alpha }`$ is the energy of the $`\alpha `$-th band of the free conduction electron problem and we have projected $`\mathrm{\Phi }_k^{}\mathrm{\Phi }_k`$ onto the eigenvectors of $`H_{kin}`$; $`\mu ,\lambda _1,\lambda _2`$ are set by saddle-point evaluations to fix the filling.

On a square lattice, C4 symmetry means that the possible overlap phases with the nearest four neighbors are $`\varphi ^\lambda =\frac{1}{2}(1,i^\lambda ,i^{2\lambda },i^{3\lambda })`$, $`\lambda =\{0,..,3\}`$. Hopping occurs in the first channel ($`\mathrm{\Gamma }=1`$) coinciding with the totally symmetric Wannier state $`\lambda =0`$. Time reversal symmetry specifies the d-wave like $`\lambda =2`$ as the second channel ($`\mathrm{\Gamma }=2`$). The form factors arise as: $`_a\varphi ^1(a)c_{i+a}^{}=\frac{1}{2}(c_{i+\widehat{x}}^{}+c_{i\widehat{y}}^{}+c_{i\widehat{x}}^{}+c_{i+\widehat{y}}^{})`$ which under fourier transform, defines $`\mathrm{\Phi }^1`$: $`\psi _{i1}^{}=_kc_k^{}e^{ikx_i}(cos(k_x)+cos(k_y))=_kc_k^{}e^{ikx_i}\mathrm{\Phi }_k^1`$. Similarly, $`\mathrm{\Phi }_k^2=c(k_y)c(k_x)`$ are found $`\mathrm{\Phi }_k^{4,3}=s(k_x)\pm is(k_y)`$ completing the basis.

The effective T$`_{c_2}`$ is found from

$`{\displaystyle \underset{k\sigma }{}}{\displaystyle \frac{f(0)f(E_{qp})}{E_{qp}}}\mathrm{\Phi }_k^2\mathrm{\Phi }_k^2+{\displaystyle \frac{2}{J}}=0`$ (14)
where $`E_{qp}=2t|b|^2\mathrm{\Phi }_k^{(1)}+\lambda _1\mu `$ is determined self-consistently subject to the conditions $`\lambda _2=0`$ and $`_kn_f(2t|b|^2\mathrm{\Phi }_k^{(1)}+\lambda _1\mu )=\frac{1|b|^2}{2}`$.

On the pyrochlore lattice, 4 bands and 6 nearest neighbors imply 24 independent symmetry channels! The requirement that the states preserve the equivalence between sites reduces this number to 6, given by

$$\varphi _\alpha ^\pm =\frac{1}{\sqrt{6}}(1,\alpha ,\alpha ^2,\pm 1,\pm \alpha ,\pm \alpha ^2)$$

where $`\alpha `$ is a cube root of unity. Of these, only the state $`\varphi ^{}=\frac{1}{\sqrt{6}}(1,1,1,1,1,1)`$ is orthogonal to the original channel and time-reversal invariant, uniquely specifying the second channel. This particular Wannier state has the property that it changes sign between neigboring tetrahedra on the lattice, as illustrated in Fig. 2. The four by four matrix $`\mathrm{\Phi }_k^{(1)}`$ now has the following parts:

$`\mathrm{\Phi }_k^{(1)}=\left(\begin{array}{cccc}0& c\left(k_a\right)& c\left(k_b\right)& c\left(k_c\right)\\ c\left(k_a\right)& 0& c\left(k_ak_b\right)& c\left(k_ak_c\right)\\ c\left(k_b\right)& c\left(k_ak_b\right)& 0& c\left(k_bk_c\right)\\ c\left(k_c\right)& c\left(k_ak_c\right)& c\left(k_bk_c\right)& 0\end{array}\right)`$ (15)
while,

$`\mathrm{\Phi }_k^{(2)}=i\left(\begin{array}{cccc}0& s\left(k_a\right)& s\left(k_b\right)& s\left(k_c\right)\\ s\left(k_a\right)& 0& s\left(k_ak_b\right)& s\left(k_ak_b\right)\\ s\left(k_b\right)& s\left(k_ak_b\right)& 0& s\left(k_bk_c\right)\\ s\left(k_c\right)& s\left(k_ak_c\right)& s\left(k_bk_c\right)& 0\end{array}\right)`$ (16)
For this case, the free hopping Hamiltonian has been treated previously by Reimers et al.berlinsky  There are two flat bands lying above two dispersing bands. Strictly speaking, since we have $`\frac{1}{2}`$ an electron per site involved and four sites, this would actually give us an insulator. However, since we have dropped the second degenerate band, a more reasonable band structure would likely include band crossings and incomplete fillings of the levels along the lines of Singh et al.singh

A large body of experimental work has been done, and we now discuss a simple interpretation of some of this work. We have four localized half-filled electron bands interacting via a Kondo-like hybridization with four itinerant levels, each less than quarter-filled, which leaves us with a hole-like Fermi surface below Tc2. In this way, our model is consistent with the observed sign change in the Hall constant which has been observed at 50 K in single crystal samples.takagi  A simple phenomenological fit consistent with the ideas we have introduced for the resistivity yields

$`\rho =\{\begin{array}{cc}\rho _i+A_hT^2\hfill & T<2K\hfill \\ \rho _o+\gamma _{ph}T+A_pT^2\hfill & T>100K\hfill \end{array}`$ (17)
where $`\rho _i=21\mu \mathrm{\Omega }cm`$ (impurity scattering in the single crystal), $`A_h=2.0\frac{\mu \mathrm{\Omega }cm}{K^2}`$ (arises from the heavy Fermi liquid) as given in ref.4, and $`\rho _o=0.27m\mathrm{\Omega }cm`$ (impurity scattering from the localized spin $`\frac{1}{2}`$), $`\gamma _{ph}=1.1\frac{\mu \mathrm{\Omega }cm}{K}`$ (phonon contribution above $`\theta _D`$) and $`A_p=2.710^3\frac{\mu \mathrm{\Omega }cm}{K^2}`$ (arising from the ’porcupine’ Fermi sea) have been extracted from Urano et al.’s datatakagi . Since $`\frac{A_h}{A_p}(\frac{m_h^{}}{m_p^{}})^2`$, this yields an effective mass ratio of $`\frac{m_h^{}}{m_p^{}}27`$. An analogous fit of the magnetic susceptibility yields

$`\chi T=\{\begin{array}{cc}\chi _hT\hfill & T<2K\hfill \\ \chi _p(T350K)+\frac{2N_A(\mu g)^2S(S+1)T}{k_B(T+\theta )}\hfill & T>400K\hfill \end{array}`$ (18)
yields $`\frac{\chi _h}{\chi _p}\frac{m_h^{}}{m_p^{}}29`$ where the $`S=\frac{1}{2}`$ ($`\theta =66K`$, g = 2.21) contribution has been first subtracted away from the high temperature susceptibility measurementskondo ; fujiwara ; takagi ; muhtar  leaving a contribution linear in T. The observed low temperature metal insulator transition is hardly surprising under pressure as broken symmetry Wannier states lie close in energy to the chosen $`\varphi ^2`$ for the second channel. Invariance of the heat capacity to 30T is indicative that the magnetic coupling J driving the spin quenching is of much higher energy than the Kondo temperature T$``$50K.

In conclusion, we have presented a physical picture within which one can understand the physics of the frustrated heavy fermion LiV2O4. A simple t-J model is capable of describing both the high temperature mixed valent phase, and the low temperature coupling between a Fermi sea and spin liquid. As no time reversal symmetry breaking has been observed, we have been able to uniquely identify the symmetry of the low temperature quasiparticles. Our picture gives a simple qualitative interpretation of the unusual resistivity, the change in sign of the Hall resistance observed and should be used to predict the result of de Haas-van Alphen experiments yet to be performed. This work should serve as a starting point for a DMFT or band structure calculation which might describe the detailed energetics involved. In particular it would be interesting to follow the compression of the lattice and see how the second channel symmetry depends on nearest neighbor distance.

We would like to acknowledge useful discussions with H. Takagi, K. LeHur and A. Rosch that led to this picture. This work was supported by National Science Foundation grant DMR 9983156.


format: text

text:
# Elliptic Flow from Minijet Production in Heavy Ion Collisions

## 1 Introduction

Ultrarelativistic heavy ion collisions provide us interesting new information on QCD under extreme conditions. Different stages of the collisions probe QCD in different regimes: early times are likely to be characterized by strong saturated gluonic fields  while creation of quark-gluon plasma (QGP) is an attractive possibility in the later stages of the collisions . To be able to understand the experimental data generated in the collisions one has to learn to disentangle between the contributions of these two (in principle) different physical processes to heavy ion observables.

The strong quasi-classical gluon fields are produced in the early stages of the collisions due to high transverse densities of color charges in the wave functions of the Lorentz-contracted nuclei . The typical color charge density is characterized by the saturation scale $`Q_s^2`$, which grows with nuclear atomic number and with energy of the collision . While being coherent over large longitudinal distances these quasi-classical fields have a rather short coherence length in the transverse direction, of the order of the inverse saturation scale $`1/Q_s`$. Saturation-inspired models have been quite successful in describing the emerging RHIC data on rapidity and centrality dependence of the total charged multiplicity of the produced particles .

In the later stages of the collisions gluons and quarks generated by the quasi-classical fields are likely to reach kinetic, and, possibly, chemical thermal equilibrium, producing the quark-gluon plasma (QGP) . To be able to study this new state of matter one should be able to distinguish the contributions to physical observables of the collective phenomena due to QGP from the effects of strong gluonic fields of the early stages. Ideally one should try to construct observables which would be sensitive only to one type of physics while being independent of the other. An example of such observables are the long range rapidity correlations predicted in , which are almost entirely due to the dynamics of initial conditions.

In this paper we would like to study the contribution of the quasi-classical gluonic fields to the elliptic flow observable $`v_2`$, defined as the second Fourier moment of the azimuthal momentum distribution of the produced particles . Elliptic flow reflects anisotropy of the transverse momentum distribution of the produced particles with respect to reaction plane. There has been a large amount of elliptic flow data produced in the heavy ion collisions at SPS  and RHIC . While the majority of these data are in good agreement with hydrodynamic simulations , the emerging new STAR data on differential elliptic flow $`v_2(p_T)`$ at high $`p_T`$ seems to deviate from hydrodynamic predictions . Instead of continuing to increase with $`p_T`$ the flow variable $`v_2(p_T)`$ saturates to approximately a constant above $`p_T=\mathrm{\hspace{0.17em}1.5}\text{GeV}`$ . The data goes up to $`p_T\mathrm{\hspace{0.17em}4.5}\text{GeV}`$ . It is natural to expect that at such high momenta the hard (perturbative) physics should be responsible for the underlying strong interactions dynamics***Very recently an attempt was made in  to explain the high-$`p_T`$ elliptic flow using pure classical gluon fields of the nuclei  yielding the elliptic flow which is too small to explain the data.. This led the authors of  to propose a model of interplay of soft and hard interactions in $`v_2`$: while the low momentum part of $`v_2`$ was still described by hydrodynamic calculations the high $`p_T`$ part was described by medium-induced radiative energy loss of partons  which would generate azimuthal anisotropy in momentum space due to coordinate space anisotropy of the overlap region. The resulting flow observable would saturate at some relatively large $`p_T`$ and then would decrease with increasing $`p_T`$.

Here we are going to propose a model of a non-flow particle production mechanism generating an elliptic flow observable $`v_2`$ which would go to a constant at large $`p_T`$ and would stay approximately constant as $`p_T`$ increases. Let us picture particle production at the early stages of a nuclear collision at high energy. In the first approximation high-$`p_T`$ particles are produced independent of each other. We can illustrate this in the framework of the saturation approach to nuclear collisions: there the gluon fields are coherent over very short transverse distances of the order of $`1/Q_s`$, with $`Q_s\mathrm{\Lambda }_{QCD}`$. The high energy wave function of a large nucleus can be viewed in the transverse plane as consisting of many independent gluon fields each of them occupying a transverse area $`1/Q_s^2`$. (In the longitudinal direction each classical field is of course coherent over the whole nucleus  and this way the nucleus can be considered as “sliced” into narrow pipes with diameter $`1/Q_s`$.) Therefore a collision with the overlap transverse area of the two nuclei $`S_{}`$ could be pictured as at least $`S_{}Q_s^2`$ independent (sub)collisions. For sufficiently central heavy ion collisions $`S_{}Q_s^2N_{part}\mathrm{\hspace{0.17em}1}`$, where $`N_{part}`$ is the number of nucleons participating in the collision. At the leading order in this large parameter $`S_{}Q_s^2`$ the two-particle multiplicity distribution factorizes into a product of two single-particle multiplicity distributions. Correlation between a pair of particles in the perturbative production mechanism happens when both particles are produced in the same subcollision, i.e., at the same impact parameter. Therefore these correlations appear at the subleading in $`S_{}Q_s^2`$ order and are suppressed by a power of $`S_{}Q_s^2`$. We are going to argue below that this does not prevent them from significantly contributing to elliptic flow observable.

The standard definition of differential elliptic flow is

$$v_2(p_T)=e^{2i(\varphi _{p_T}\mathrm{\Phi }_R)}=\mathrm{cos}2(\varphi _{p_T}\mathrm{\Phi }_R)$$
(1)
where $`\varphi _{p_T}`$ is the azimuthal angle of the produced particle with the value of transverse momentum $`p_T`$, $`\mathrm{\Phi }_R`$ is the azimuthal angle of the reaction plane and the brackets denote statistical averaging over different events which leaves non-zero only the contribution of the cosine in Eq. (1). The reaction plane is determined by averaging over particles produced in a heavy ion collision with certain weights designed to optimize the reaction plane resolution . Let us imagine that the particle with momentum $`p_T`$ in Eq. (1) was produced in the same subcollision with some other particle, which contributed to determination of the reaction plane angle $`\mathrm{\Phi }_R`$. Then the particle $`p_T`$ would be correlated with the reaction plane by these non-flow azimuthal correlations. The averaged over events correlations do not disappear at high $`p_T`$, as will be shown below. Therefore the elliptic flow observable $`v_2`$ may be potentially sensitive to the non-flow correlations originating from minijet production in the initial conditions. To avoid this problem the authors of  (see also ) introduced a cumulant approach to flow analysis. In terms of the saturation model contribution of minijets to higher order cumulants defined in  is suppressed by powers of $`S_{}Q_s^2N_{part}`$. This led the authors of  to suggest that the high-cumulant flow analysis would be relatively free of minijet effects compared to standard flow analysis.

There are certain potential dangers in this conclusion which have to be addressed in an explicit calculation. One is that for peripheral collisions $`S_{}Q_s^2N_{part}`$ is not such a large number and therefore minijet effects should not be suppressed anymore. From the experimental data  we know that elliptic flow is strongest in peripheral collisions making the minijet contribution also large. Another potentially dangerous question is whether the parameter $`S_{}Q_s^2`$ appears in actual calculations with some numerically small prefactor making inverse powers of $`S_{}Q_s^2`$ not too small. After all the observed differential $`v_2`$ is also not a very large quantity being of the order of $`1015\%`$ at RHIC and may be very sensitive to such corrections.

To clarify the questions mentioned above we are going to perform an explicit calculation of the minijet contribution to the elliptic flow observable. The paper is structured as follows. In Sect. 2 we show that the standard determination of flow employing Eq. (10) and the flow extracted from the 2nd cumulant proposed in  should yield the same result for $`v_2`$ for a trivial choice of weights. This result was confirmed by an explicit flow analysis using two different methods at STAR giving almost the same result for $`v_2`$ .

The second order cumulant from  is nothing but a two-particle correlation function, a contribution to which from the minijet production is calculated in Sect. 3. To construct a model of minijet production which incorporates both hard and soft physics we have used $`k_T`$ factorized expression for two-jet production with the unintegrated gluon distributions given by the quasi-classical Glauber-Mueller expression . These gluon distributions are characterized by the saturation scale $`Q_s\mathrm{\Lambda }_{QCD}`$ which insures applicability of small coupling approaches down to rather small transverse momenta of the produced particles. Our model of course does not yield us the exact two-particle production cross section, but gives a realistic approximation similar to the one used in  to describe the multiplicity distributions at RHIC. An exact calculation of double inclusive minijet production in the quasi-classical framework appears to be rather complicated and is left for further research: even the single inclusive cross section has not been yet unambiguously theoretically determined, despite the extensive efforts .

In Sect. 3 we calculate the differential flow observable resulting from these two-particle correlations. Our final result is given in Eq. (40). The obtained differential elliptic flow $`v_2(p_T)`$ starts increasing as a power of $`p_T`$ for small $`p_TQ_s`$ (see Eq. (41)) and then saturates to a slow logarithmic growth for $`p_TQ_s`$ (see Eq. (42)). We also derive the centrality dependence of the minijet flow contribution (see Eq. (45)).

In Sect. 4 after making some simple assumptions about the gluon distributions employed in the minijet production cross section we fit the differential elliptic flow STAR data  using the flow from Eq. (40) with $`\alpha _s=\mathrm{\hspace{0.17em}0.3}`$, $`Q_s=\mathrm{\hspace{0.17em}1}\text{GeV}^2`$. These values are in agreement with the ones used in the saturation-inspired analysis of the multiplicity data in . The fit is shown in Fig. 3. We can also fit the centrality dependence of $`v_2`$ with our minijet model by noting that the integrated flow should scale as $`v_2(B)\mathrm{\hspace{0.17em}1}/\sqrt{S_{}Q_s^2}\mathrm{\hspace{0.17em}1}/\sqrt{N_{part}}`$ (see Fig. 4).

We discuss the results in Sect. 5 by stating that while a more involved numerical analysis is still needed to analyze the emerging RHIC data on elliptic flow we have demonstrated that the contribution of minijets to the standard flow analysis is very large, possibly accounting for most of high-$`p_T`$ data. Thus it appears that the standard flow analysis is heavily “contaminated” by minijets which prevent direct measurements of the contribution of collective QGP effects to elliptic flow. At the same time flow analysis seems to be rather sensitive to details of saturation physics and could be used for determination of nuclear saturation scales.

## 2 Different Methods of Flow Analysis

In  Borghini et al proposed a new and interesting approach to flow analysis employing higher order cumulants. Let us outline some important features of the approach presented in  for the lowest order cumulant. The second order cumulant is defined for two particles with azimuthal angles $`\varphi _1`$ and $`\varphi _2`$ as

$$e^{2i(\varphi _1\varphi _2)}=e^{2i(\varphi _1\varphi _2)}e^{2i\varphi _1}e^{2i\varphi _2}.$$
(2)
The average $`e^{2i\varphi _1}`$ vanishes since the angle here is measured in the laboratory and heavy ion collisions are azimuthally symmetric after averaging over many events. For flow correlations this means that the angle in $`e^{2i\varphi _1}`$ is not measured with respect to reaction plane, but with respect to some fixed direction in the detector. Assuming that the particles $`1`$ and $`2`$ are correlated with each other only through the flow correlations with the reaction plane the authors of  wrote

$$e^{2i(\varphi _1\varphi _2)}=e^{2i(\varphi _1\varphi _2)}=e^{2i(\varphi _1\mathrm{\Phi }_R)}e^{2i(\mathrm{\Phi }_R\varphi _2)}=(v_2)^2,$$
(3)
where the definition of elliptic flow from Eq. (1) was employed. A new method for measuring elliptic flow was proposed in  using the two particle (and higher order) cumulants of Eq. (3). If we fix the transverse momentum of particle $`1`$ to be $`p_T`$ and average over all momenta of the particle $`2`$ over many events then as one can see from Eq. (3)

$$e^{2i(\varphi _1(p_T)\varphi _2)}=v_2(p_T)v_2$$
(4)
with $`v_2`$ the elliptic flow variable averaged over all $`p_T`$. At the same time if we do not impose any restrictions on the transverse momenta of both particles we get

$$e^{2i(\varphi _1\varphi _2)}=v_2^2.$$
(5)
From Eqs. (4) and (5) noting that only cosine components of the exponents survive the averaging we obtain the following expression for differential elliptic flow

$$v_2(p_T)=\frac{\mathrm{cos}(2(\varphi _1(p_T)\varphi _2))}{\sqrt{\mathrm{cos}(2(\varphi _1\varphi _2))}}.$$
(6)
Let us define the event-averaged two-particle multiplicity distribution function

$$P(k_1,k_2,\underset{¯}{B})=\frac{dN}{d^2k_1dy_1d^2k_2dy_2}(\underset{¯}{B}),$$
(7)
where the transverse momenta of the particles are $`\underset{¯}{k}_1`$ and $`\underset{¯}{k}_2`$, while $`y_1`$ and $`y_2`$ are their rapidities and we average over all events with the impact parameter $`\underset{¯}{B}`$ between the two nuclei. This distribution can be written as a sum of uncorrelated and correlated terms

$$P(k_1,k_2,\underset{¯}{B})=\frac{dN}{d^2k_1dy_1}\frac{dN}{d^2k_2dy_2}+\frac{dN_{corr}}{d^2k_1dy_1d^2k_2dy_2},$$
(8)
where we suppressed the impact parameter dependence. The correlated term in Eq. (8) is usually much smaller than the uncorrelated one as it is suppressed by a power of $`Q_s^2S_{}N_{part}`$. Using Eq. (7) we rewrite Eq. (6) as

$`v_2(k_1,\underset{¯}{B})={\displaystyle \frac{d^2k_2𝑑y_2𝑑\varphi _1𝑑y_1P(k_1,k_2,\underset{¯}{B})\mathrm{cos}(2(\varphi _1\varphi _2))}{d^2k_2𝑑y_2𝑑\varphi _1𝑑y_1P(k_1,k_2,\underset{¯}{B})}}`$
$$\times \left(\frac{d^2k_1𝑑y_1d^2k_2𝑑y_2P(k_1k_2,\underset{¯}{B})}{d^2k_1𝑑y_1d^2k_2𝑑y_2P(k_1,k_2,\underset{¯}{B})\mathrm{cos}(2(\varphi _1\varphi _2))}\right)^{1/2},$$
(9)
where we have relabeled the transverse momenta of the particles to be $`\underset{¯}{k}_1`$ and $`\underset{¯}{k}_2`$. Eq. (9) gives us a way of calculating elliptic flow from two-particle correlation functions. Higher order cumulants would yield us ways of calculating flow from higher order particle correlations , though one can not go to arbitrary high order cumulants due to lack of statistics there. An analysis of RHIC data has been performed at STAR  using both conventional and cumulant approaches explicitly demonstrating that the flow obtained by the second cumulant technique of Eq. (9) is consistent with the flow extracted using conventional techniques. That is the approach of Eq. (9) is equivalent to the conventional flow analysis. To see this let us first note that in the actual standard flow analysis one has to take into account the resolution in the event plane determination . This would modify Eq. (1) giving

$$v_2(p_T)=\frac{\mathrm{cos}(2(\varphi _{p_T}\mathrm{\Psi }_R))}{\sqrt{2\mathrm{cos}(2(\mathrm{\Psi }_R^a\mathrm{\Psi }_R^b))}}$$
(10)
where $`\mathrm{\Psi }_R^a`$ and $`\mathrm{\Psi }_R^b`$ are the event plane angles determined in two different sub-events labeled $`a`$ and $`b`$ with the particle multiplicity $`N/2`$ in each of them while $`\mathrm{\Psi }_R`$ is the event plane angle determined in the full event with multiplicity $`N`$ . Eq. (10) is quite similar to Eq. (6). To show that the two equations are almost equivalent let us go back to exponential notation and rewrite the numerator of Eq. (6) as

$$\mathrm{cos}(2(\varphi _1(p_T)\varphi _2))=e^{2i(\varphi _1(p_T)\varphi _2(k_T))}_{allk_Tp_T}_{events}$$
(11)
where we first average over all particles (with various $`k_T`$) in a given event other than the chosen particle with momentum $`p_T`$, just like in the standard flow analysis , and then average over all events. On the other hand, the azimuthal angle of the reaction plane $`\mathrm{\Psi }_R`$ in Eq. (10) can be defined by

$$Q_Ne^{2i\mathrm{\Psi }_R}=e^{2i\varphi _{k_T}}_{allk_Tp_T},$$
(12)
where $`Q_N`$ is a real number. If one neglects correlations between the particles and multiplies Eq. (12) by its complex conjugate avergaing over events with the same total multiplicity $`N`$ one easily gets $`<Q_N^2>=\mathrm{\hspace{0.17em}1}/N`$ . Therefore at the leading order in $`N`$ we have (for large $`N`$) $`Q_N\mathrm{\hspace{0.17em}1}/\sqrt{N}`$. In Eq. (12) we for the moment forget about the transverse momentum cutoffs imposed in event plane determination and subtleties related to different choices of weights and detector imperfections. Substituting Eq. (12) into the numerator of Eq. (10) we obtain

$`\mathrm{cos}(2(\varphi _{p_T}\mathrm{\Psi }_R))_{events}=e^{2i(\varphi _{p_T}\mathrm{\Psi }_R)}_{events}={\displaystyle \frac{1}{Q_N}}e^{2i(\varphi _{p_T}\varphi _{k_T})}_{allk_Tp_T}_{events}`$
$$\frac{1}{Q_N}e^{2i(\varphi _{p_T}\varphi _{k_T})}_{allk_Tp_T}_{events}$$
(13)
where we assumed that for a fixed impact parameter collisions the particle multiplicity $`N`$ (and, therefore, $`Q_N`$) is independent of two-particle correlations. Eq. (13) is identical to Eq. (11) up to a factor of $`1/Q_N`$. Similarly one can show that the denominators in Eqs. (10) and (6) are equal to each other up to a factor of $`\sqrt{2/Q_{N/2}^2}`$:

$`2\mathrm{cos}(2(\mathrm{\Psi }_R^a\mathrm{\Psi }_R^b))_{events}=\mathrm{\hspace{0.17em}2}{\displaystyle \frac{1}{Q_{N/2}^2}}e^{2i(\varphi _{k_{T1}}\varphi _{k_{T2}})}_{alldifferentk_{T1},k_{T2}}_{events}`$
$$\frac{2}{Q_{N/2}^2}\mathrm{cos}(2(\varphi _1\varphi _2)).$$
(14)
Recalling that for large multiplicities $`Q_N1/\sqrt{N}`$ we see that at the leading order in $`N`$: $`\sqrt{Q_{N/2}^2}/\sqrt{2}Q_NQ_{N/2}/\sqrt{2}Q_N1`$. Extra factors in the numerator and denominator of Eq. (10) cancel reducing it to Eq. (6). Thus we showed that Eqs. (10) and (6) are identical in the limit of large multiplicity $`N`$. We have proven that the standard definition of flow from Eq. (10) is equivalent to the definition introduced in  given here by Eq. (6)). Therefore in order to calculate the differential elliptic flow $`v_2(p_T)`$ we need only to calculate the two-particle correlation function $`P(k_1,k_2,\underset{¯}{B})`$ and substitute it into Eq. (9).

## 3 Calculation of $`v_2`$

To calculate the two-particle correlation coefficient in Eq. (8) one needs to know single and double particle multiplicity distributions. While the quasi-classical single particle production mechanisms in heavy ion collisions have been extensively studied , the double inclusive particle production has not been calculated yet. One needs to calculate production of a pair of particles with comparable transverse momenta $`|\underset{¯}{k}_1||\underset{¯}{k}_2|Q_s`$ and rapidities $`y_1y_2`$ in the framework of McLerran-Venugopalan model , that is resumming all powers of $`Q_s^2/k_{}^2`$ . Unlike the single gluon production this process involves an extra gluon and can not be described by the classical field methods of . At the lowest order in $`Q_s^2/k_{}^2`$ the two gluon production amplitude is equivalent to the real part of NLO BFKL kernel  and is rather sophisticated . Going beyond leading order in $`Q_s^2/k_{}^2`$ appears to be tremendously complicated and we will not address this problem here. Instead we are going to construct a model of correlated two-particle production employing quasi-classical gluon distributions from  in the production formula inspired by collinear factorization , similar to how it was done in  in describing the multiplicity data at RHIC. We would also assume that in any given event either $`y_1y_2`$ or $`y_1y_2`$. Of course this assumption does not hold in actual flow analyses  and we are making it just to simplify the calculations. While, as discussed above, a more detailed calculation would still be required to obtain an exact expression for the correlation coefficient in Eq. (8), we believe that it would only introduce numerical corrections to our approach leaving qualitative results the same.

Consider inclusive production of one or two gluons in scattering of two quarks at high energy as shown in Fig. 1, where the blobs denote effective Lipatov vertices. It is convenient to use the Sudakov decomposition of gluons’ four-momenta:

$$q_i=\alpha _ip_1+\beta _ip_2+\underset{¯}{q}_i,i=1,2,3.$$
(15)
In Fig. 1(a) the dominant contribution stems from the multi-regge kinematical region where $`\alpha _1\alpha _2`$ and $`\beta _2\beta _1`$. Using Ward identities the $`t`$-channel gluon propagators can be written as

$$D^{\mu \nu }(q_i)=\frac{2}{s}\frac{\underset{¯}{q}_i^\mu \underset{¯}{q}_i^\nu }{\alpha _i\beta _i},$$
(16)
where $`s=\mathrm{\hspace{0.17em}2}(p_1p_2)`$ is center-of-mass energy squared. With this representation of the $`t`$-channel gluon propagator effective Lipatov’s vertex  for $`s`$-channel gluon production (black blobs in the Fig. 1) reduces to the usual three-gluon vertex $`\mathrm{\Gamma }_{\mu \nu }^\rho `$ and can be easily calculated:

$$L^\rho =\underset{¯}{q}_1^\mu \underset{¯}{q}_2^\nu \mathrm{\Gamma }_{\mu \nu }^\rho =\frac{1}{2}\alpha _1\beta _2s\left[\left(\alpha _1+\frac{2\underset{¯}{q}_1^2}{\beta _2s}\right)p_1^\rho +\left(\beta _2+\frac{2\underset{¯}{q}_2^2}{\alpha _1s}\right)p_2^\rho (\underset{¯}{q}_1+\underset{¯}{q}_2)^\rho \right].$$
(17)
In our calculation we will need the square of Lipatov’s vertex given by

$$(L^\rho )^2=\alpha _1\beta _2s\underset{¯}{q}_1^2\underset{¯}{q}_2^2.$$
(18)
Using Eq. (16), Eq. (17), Eq. (18) and the four-dimensional momentum element decomposition

$$d^4q_i=\frac{s}{2}d^2\underset{¯}{q}_id\alpha _id\beta _i$$
(19)
we evaluate the single and double distribution of the produced particles shown in Fig. 1

$$\frac{dN}{d^2k_1dy_1}=\frac{\alpha _s^3}{\pi ^2}\frac{2C_FA^2}{S_{}\underset{¯}{k}_1^2}d^2q_1\frac{1}{\underset{¯}{q}_1^2(\underset{¯}{k}_1\underset{¯}{q}_1)^2},$$
(20)
$$\frac{dN_{corr}}{d^2k_1dy_1d^2k_2dy_2}=\frac{\alpha _s^4}{\pi ^4}\frac{N_cC_FA^2}{S_{}\underset{¯}{k}_1^2\underset{¯}{k}_2^2}d^2q_1\frac{1}{\underset{¯}{q}_1^2(\underset{¯}{k}_1+\underset{¯}{k}_2\underset{¯}{q}_1)^2},$$
(21)
where we changed variables to $`\underset{¯}{k}_1=\underset{¯}{q}_1\underset{¯}{q}_2`$, $`\underset{¯}{k}_2=\underset{¯}{q}_2\underset{¯}{q}_3`$ (the momenta of produced gluons). Throughout the paper we assume for simplicity that the nucleus has a cylindrical shape with radius $`R`$ and “height” $`2R`$ and the axis of the cylinder coincides with the collision axis. $`S_{}=S_{}(\underset{¯}{B})`$ is the nucleus overlap transverse area which depends on nuclei relative impact parameter $`\underset{¯}{B}`$. We assume that colliding nuclei are identical having atomic number $`A`$ each.

In fact colliding quarks in Fig. 1 are not free but confined inside nucleons. If momenta of produced gluons are higher than all scales characterizing the nucleus (including $`Q_s`$) we can use $`k_t`$-factorization to write the distributions from Eq. (20) and Eq. (21) as convolutions of hard amplitudes with parton distributions in nuclei which are predominantly gluonic in the high energy region $`st`$. Employing the relation

$$xG_A(x,\underset{¯}{q}^2)=AxG(x,\underset{¯}{q}^2)=A\frac{\alpha _sC_F}{\pi }\mathrm{ln}(\underset{¯}{q}^2/\mu ^2),$$
(22)
with $`\mu `$ some non-perturbative cutoff the one and two-particle distributions can be written as

$$\frac{dN}{d^2k_1dy_1}=\frac{2\alpha _s}{C_FS_{}}\frac{1}{\underset{¯}{k}_1^2}d^2q_1\frac{dxG_A}{d\underset{¯}{q}_1^2}\frac{dxG_A}{d(\underset{¯}{k}_1\underset{¯}{q}_1)^2}$$
(23)
$$\frac{dN_{corr}}{d^2k_1dy_1d^2k_2dy_2}=\frac{N_c\alpha _s^2}{\pi ^2C_FS_{}}\frac{1}{\underset{¯}{k}_1^2\underset{¯}{k}_2^2}d^2q_1\frac{dxG_A}{d\underset{¯}{q}_1^2}\frac{dxG_A}{d(\underset{¯}{k}_1+\underset{¯}{k}_2\underset{¯}{q}_1)^2}.$$
(24)
As was argued in the Introduction, if the number of partons in each nuclei is large enough then their mutual interactions must be taken into account. This may happen either due to increase of parton density in each nucleon by subsequent emission of gluons in course of quantum evolution as collision energy increases  or due to enhancement of nucleon parton density in a large nucleus by the atomic number $`A`$ . The latter effect leads to creation of the non-Abelian Weiszäcker-Williams (WW) field $`\underset{¯}{A}^{WW}(\underset{¯}{z})`$ of the nucleus which is a strong classical gluon field already at moderate energies . The WW field gives rise to the unintegrated gluon distribution given by

$`{\displaystyle \frac{dxG_A(x,\underset{¯}{q}^2)}{d\underset{¯}{q}^2}}`$ $`=`$ $`{\displaystyle \frac{2}{(2\pi )^2}}{\displaystyle d^2\underset{¯}{z}e^{i\underset{¯}{z}\underset{¯}{q}}d^2\underset{¯}{b}\mathrm{Tr}\underset{¯}{A}^{WW}(\underset{¯}{0})\underset{¯}{A}^{WW}(\underset{¯}{z})}`$ (25)
$`=`$ $`{\displaystyle \frac{2}{\pi (2\pi )^2}}{\displaystyle d^2\underset{¯}{z}e^{i\underset{¯}{z}\underset{¯}{q}}\frac{S_{}C_F}{\alpha _s\underset{¯}{z}^2}\left(1e^{\frac{1}{4}\underset{¯}{z}^2Q_s^2}\right)},`$
where $`\underset{¯}{b}`$ is the gluon’s impact parameter (which we can trivially integrate over in a cylindrical nucleus) and

$$Q_s^2(\underset{¯}{z})=\frac{4\pi ^2\alpha _sN_c}{N_c^21}\rho xG(x,1/\underset{¯}{z}^2)T(\underset{¯}{b}),$$
(26)
with $`\rho =A/[2\pi R^3]`$ the atomic number density in the cylindrical nucleus with atomic number $`A`$, and $`T(\underset{¯}{b})`$ the nucleus profile function equal to $`2R`$ for the cylindrical nucleus considered here. This provides the initial condition to the nonlinear quantum evolution of the gluon distribution with energy in the high parton density region . $`Q_s`$ is a scale at which nonlinear nature of the gluon field becomes evident. We suggest using the classical expression Eq. (25) as an approximation to the exact gluon field of the nucleus. This is a justified approximation as long as $`\alpha _s\mathrm{ln}s\genfrac{}{}{0pt}{}{_<}{}\mathrm{\hspace{0.17em}1}`$, i.e., when corrections due to quantum evolution are small.

It is usually assumed that $`k_T`$-factorization holds in high parton density regime as well as in the linear one . Phenomenological models for heavy ion collisions which employ this assumption together with nonlinear evolution for gluon distributions proved to be successful in describing experimental data at SPS and RHIC . This implies that the use of $`k_T`$-factorization is a quite good approximation for the SPS and RHIC kinematical regions even for $`k_T<Q_s`$. Therefore substituting Eq. (25) into Eq. (23) and into Eq. (24) we obtain the following expressions for the single and double gluon distributions

$$\frac{dN}{d^2k_1dy_1}=\frac{C_FS_{}}{\alpha _s\underset{¯}{k}_1^2}\frac{4K_1}{\pi ^3}_0^{\mathrm{}}\frac{dz}{z^3}J_0(k_1z)\left(1e^{\underset{¯}{z}^2Q_s^2/4}\right)^2,$$
(27)
$$\frac{dN_{corr}}{d^2k_1dy_1d^2k_2dy_2}=\frac{N_cC_FS_{}}{\underset{¯}{k}_1^2\underset{¯}{k}_2^2}\frac{K_2}{\pi ^6}d^2\underset{¯}{z}\frac{1}{\underset{¯}{z}^4}e^{i\underset{¯}{z}(\underset{¯}{k}_1+\underset{¯}{k}_2)}\left(1e^{\underset{¯}{z}^2Q_s^2/4}\right)^2,$$
(28)
where the saturation scales in both nuclei are the same since the nuclei are cylindrical and identical. Generalization of Eqs. (27) and (28) to a spherical nucleus is straightforward.

We have to point out again that the single particle distribution is, in principle, known better than displayed here in Eq. (27) . For instance the distribution in Eq. (27) is not infrared safe, while the correct distribution derived in  has no infrared divergences. However, as we need both single and double particle distributions to obtain elliptic flow using Eq. (9) we should calculate both of them in the framework of the same model. Since the exact calculation of the double gluon distribution does not seem feasible at the moment we have to calculate both distributions in the same $`k_T`$-factorization approach inserting $`K`$-factors $`K_1`$ and $`K_2`$ to correct the normalization of the approximation to include higher order corrections , which we have done in Eqs. (27) and (28). The value of the $`K`$-factors will be fixed later. We will determine $`K_1`$ by comparing particle multiplicity per unit rapidity ($`dN/dy`$) resulting from Eq. (27) to the total multiplicity observed at RHIC . To fix $`K_2`$ we consider two-particle production cross section, which is proportional to the two-particle multiplicity distribution function $`P(k_1,k_2,\underset{¯}{B})`$ from Eq. (8). In the limit of large transverse momentum $`k_1k_2p_T`$ the second term in Eq. (8) given by Eq. (28) falls off at most as $`1/p_T^6`$ and has a collinear singularity $`\underset{¯}{k}_1+\underset{¯}{k}_2=0`$. Therefore it dominates over the first term in Eq. (8) given by Eq. (27) squared, which gives a $`1/p_T^8`$ fall off of the cross section. For large $`\underset{¯}{k}_1=\underset{¯}{k}_2`$ Eq. (8) and, therefore, Eq. (28) should match onto the corresponding collinear factorization expression for back-to-back jets . To fix the normalization we expand the term in the parentheses of Eq. (28) to the lowest order and integrate over one of the transverse momenta obtaining

$$\frac{d\sigma _{corr}}{dp_T^2dy_1dy_2}=K_2\frac{9\pi \alpha _s^2}{4p_T^4}[xG(x,p_T^2)]^2$$
(29)
in agreement with collinear factorization result at mid-rapidity . We put $`N_c=3`$ explicitly in Eq. (29). Collinear factorization models  are rather successful in describing the high-$`p_T`$ particle spectra in hadronic and heavy ion collisions when putting $`K_2=2`$ to correct the lowest order perturbative expression for next-to-leading order effects. For our model to be in agreement with these high-$`p_T`$ data we will have to also put $`K_2=2`$ when trying to describe the flow data in the next section. We have to note that since the gluon distributions in Eq. (25) do not include DGLAP  evolution in them and Eq. (28) does not have jet quenching effects  in it our model can not be applied to describe the high-$`p_T`$ spectra and a complete DGLAP-evolved gluon distributions along with jet quenching should be used in a more detailed numerical treatment of the problem as was done for particle spectra in .

Substituting Eqs. (27) and (28) into Eq. (8) we are ready to calculate $`v_2(p_T,\underset{¯}{B})`$. To calculate the integral in the numerator of the first line of Eq. (9) we note that only the correlated part of $`P(k_1,k_2,\underset{¯}{B})`$ (the second term in Eq. (8)) contributes there. Let us denote the angle between two-vectors $`\underset{¯}{z}`$ and $`\underset{¯}{k}_1`$ by $`\theta `$ and the angle between $`\underset{¯}{k}_2`$ and $`\underset{¯}{k}_1`$ by $`\varphi _2`$. The third angle defining direction of $`\underset{¯}{k}_1`$ is free, so we can just integrate over it obtaining a factor of $`2\pi `$. Plugging Eq. (28) into Eq. (8) and using the latter in Eq. (9) we get

$`{\displaystyle d^2k_2𝑑\varphi _1P(k_1,k_2,\underset{¯}{B})\mathrm{cos}(2(\varphi _1\varphi _2))}={\displaystyle d^2k_2𝑑\varphi _1\frac{dN_{corr}}{d^2k_1dy_1d^2k_2dy_2}\mathrm{cos}(2(\varphi _1\varphi _2))}`$
$`={\displaystyle \frac{N_cC_FS_{}}{\underset{¯}{k}_1^2}}{\displaystyle \frac{K_2}{\pi ^5}}{\displaystyle \frac{dz}{z^3}𝑑\theta \frac{dk_2^2}{k_2^2}\left(1e^{\underset{¯}{z}^2Q_s^2/4}\right)^2_0^{2\pi }𝑑\varphi _2e^{izk_2\mathrm{cos}(\varphi _2\theta )}e^{izk_1\mathrm{cos}\theta }\mathrm{cos}2\varphi _2}`$
$`={\displaystyle \frac{N_cC_FS_{}}{\underset{¯}{k}_1^2}}{\displaystyle \frac{4K_2}{\pi ^3}}{\displaystyle \frac{dz}{z^3}\frac{dk_2^2}{k_2^2}\left(1e^{\underset{¯}{z}^2Q_s^2/4}\right)^2J_2(k_1z)J_2(k_2z)}`$
$`={\displaystyle \frac{N_cC_FS_{}}{\underset{¯}{k}_1^2}}{\displaystyle \frac{4K_2}{\pi ^3}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dz}{z^3}}\left(1e^{\underset{¯}{z}^2Q_s^2/4}\right)^2J_2(k_1z).`$ (30)
To calculate the integral in the denominator of the first line of Eq. (9) we first note that the final state multiplicity per unit rapidity in heavy ion collisions was calculated in  (see also ) and reads

$$\frac{dN}{dy}=c\frac{S_{}C_FQ_s^2}{\alpha _s\mathrm{\hspace{0.17em}2}\pi ^2}$$
(31)
with $`c`$ the gluon liberation coefficient which was calculated in  to be $`c=2\mathrm{ln}2`$. Following the same steps as in derivation of Eq. (30) we find contributions to the integral in the denominator of the first line of Eq. (9) of uncorrelated two-particle distribution

$`{\displaystyle d^2k_2𝑑\phi _1\frac{dN}{d^2k_1dy_1}\frac{dN}{d^2k_2dy_2}}={\displaystyle \frac{dN}{dy_2}}{\displaystyle \frac{dN}{k_1dk_1dy_1}}`$
$`=c{\displaystyle \frac{C_F^2S_{}^2Q_s^2}{\alpha _s^2\underset{¯}{k}_1^2}}{\displaystyle \frac{4K_1}{\pi ^4}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dz}{z^3}}J_0(k_1z)\left(1e^{\underset{¯}{z}^2Q_s^2/4}\right)^2`$ (32)
and of correlated two-particle distribution

$`{\displaystyle d^2k_2𝑑\phi _1\frac{dN_{corr}}{d^2k_1dy_1d^2k_2dy_2}}`$
$`={\displaystyle \frac{N_cC_FS_{}}{\underset{¯}{k}_1^2}}{\displaystyle \frac{8K_2}{\pi ^3}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dz}{z^3}}\left(1e^{\underset{¯}{z}^2Q_s^2/4}\right)^2J_0(k_1z)\mathrm{ln}{\displaystyle \frac{1}{z\mu }}.`$ (33)
As one can obviously see the integral in Eq. (32) is enhanced by factor of $`S_{}Q_s^2/\alpha _s^2`$ as compared to the integral in Eq. (33). This means that the total number of correlated pairs of particles is negligible compared to the number of the uncorrelated pairs given by the square of the total multiplicity. Therefore the correlated contribution form Eq. (33) can be neglected compared to the uncorrelated contribution in Eq. (32).

To estimate the integrals in the second line of Eq. (9) we need to integrate Eq. (30) over momenta $`k_1`$

$`{\displaystyle d^2k_1d^2k_2P(k_1,k_2,\underset{¯}{B})\mathrm{cos}(2(\varphi _1\varphi _2))}`$
$`={\displaystyle _0^{\mathrm{}}}𝑑k_1k_1{\displaystyle \frac{N_cC_FS_{}}{k_1^2}}{\displaystyle \frac{4K_2}{\pi ^3}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{dz}{z^3}}\left(1e^{\underset{¯}{z}^2Q_s^2/4}\right)^2J_2(k_1z){\displaystyle \frac{N_cC_FS_{}Q_s^2K_2\mathrm{ln}2}{2\pi ^3}},`$ (34)
where we assumed that $`Q_s^2`$ of Eq. (26) is approximately independent of transverse coordinate $`\underset{¯}{z}`$ neglecting the logarithm of Eq. (22). To complete the calculation of the integrals in Eq. (9) we note that

$$d^2k_1d^2k_2P(k_1,k_2,\underset{¯}{B})=\frac{dN}{dy_1}\frac{dN}{dy_2}=\left(c\frac{S_{}C_FQ_s^2}{\alpha _s\mathrm{\hspace{0.17em}2}\pi ^2}\right)^2.$$
(35)
Inserting Eqs. (30), (32), (3) and (35) into Eq. (9) and noting that the integration over rapidities is trivial and cancels out between different integrals we obtain

$$v_2(p_T,\underset{¯}{B})=\alpha _s\left(\frac{\pi N_cK_2}{2\mathrm{ln}2C_FS_{}Q_s^2K_1^2}\right)^{1/2}\frac{_0^{\mathrm{}}\frac{dz}{z^3}J_2(p_Tz)\left(1e^{\underset{¯}{z}^2Q_s^2/4}\right)^2}{_0^{\mathrm{}}\frac{dz}{z^3}J_0(p_Tz)\left(1e^{\underset{¯}{z}^2Q_s^2/4}\right)^2}.$$
(36)
We have changed the transverse momentum back to $`p_T`$ to comply with conventional notation. Eq. (36) gives the minijet contribution to the differential elliptic flow for collision of two identical nuclei at given impact parameter $`\underset{¯}{B}`$. $`v_2(p_T)`$ measured in experiments is averaged over all impact parameters

$$v_2(p_T)=\frac{d^2Bv_2(p_T,\underset{¯}{B})(dN/d^2p_Tdy)}{d^2B(dN/d^2p_Tdy)}.$$
(37)
In our model each nucleus is cylindrical, so that the dependence on the impact parameter $`B`$ comes only from the nuclear overlap area $`S_{}(B)`$. The overlap area is

$$S_{}=R^2(\beta \mathrm{sin}\beta )=R^2(2\mathrm{arccos}(B/2R)\mathrm{sin}(2\mathrm{arccos}(B/2R)),$$
(38)
where $`\beta `$ is the opening angle in the transverse plane as shown in Fig. 2.

From Eq. (27) we see that $`dN/d^2p_TdyS_{}(\underset{¯}{B})`$ and from Eq. (36) one concludes that $`v_2(p_T,\underset{¯}{B})\mathrm{\hspace{0.17em}1}/\sqrt{S_{}(\underset{¯}{B})}`$. Therefore after using Eqs. (27) and (36) in Eq. (37) the impact parameter averaging reduces to

$$\frac{_0^{2R}𝑑BB\sqrt{2\mathrm{arccos}(B/2R)\mathrm{sin}(2\mathrm{arccos}(B/2R)}}{R_0^{2R}dBB(2\mathrm{arccos}(B/2R)\mathrm{sin}(2\mathrm{arccos}(B/2R))}\frac{.99}{R}\sqrt{\frac{\pi }{S_{}^A}},$$
(39)
where $`S_{}^A=\pi R^2`$ is the cross sectional area of the nucleus. The final result for $`v_2(p_T)`$ reads

$$v_2(p_T)=\alpha _s\left(\frac{\pi ^2N_cK_2}{2\mathrm{ln}2C_FS_{}^AQ_s^2K_1^2}\right)^{1/2}\frac{_0^{\mathrm{}}\frac{dz}{z^3}J_2(p_Tz)\left(1e^{\underset{¯}{z}^2Q_s^2/4}\right)^2}{_0^{\mathrm{}}\frac{dz}{z^3}J_0(p_Tz)\left(1e^{\underset{¯}{z}^2Q_s^2/4}\right)^2}.$$
(40)
Eq. (40) is our main result for differential elliptic flow. Let us first study the qualitative features of the flow from Eq. (40) before using it to fit RHIC data in Sect. 4.

In the small momentum region, $`p_TQ_s`$, we can neglect the exponents in the numerator and denominator of Eq. (40). Expanding the Bessel functions for small $`p_Tz`$ and cutting off the integrals over $`z`$ by $`1/p_T`$ from above and $`1/Q_s`$ from below we end up with

$$v_2(p_T)\alpha _s\left(\frac{\pi ^2N_cK_2}{2\mathrm{ln}2C_FS_{}^AQ_s^2K_1^2}\right)^{1/2}\frac{p_T^2}{4Q_s^2}\mathrm{ln}\frac{Q_s}{p_T},p_TQ_s.$$
(41)
We see that minijet $`v_2`$ is an increasing function of transverse momentum, which is in qualitative agreement with the data.

The asymptotic behavior of $`v_2`$ at high transverse momenta of the particles, $`p_TQ_s`$, can be found by expanding the exponents in the numerator and denominator of Eq. (40). Here it is essential to keep logarithms arising from Eqs. (26) and (22). Simple integration yields

$$v_2(p_T)\alpha _s\left(\frac{\pi ^2N_cK_2}{2\mathrm{ln}2C_FS_{}^AQ_s^2K_1^2}\right)^{1/2}\mathrm{ln}(p_T/\mu ),p_TQ_s.$$
(42)
Therefore we conclude that while the differential elliptic flow given by Eq. (40) increases for small transverse momenta, at $`p_TQ_s`$ it turns over and its growth slows down to just a logarithmic increase. As the transverse momentum increases the DGLAP  logarithms should become important modifying Eq. (42). The logarithm in Eq. (42) might also be an artifact of our approach and it might disappear if a more detailed analysis is carried out with inclusion of evolution effects in the gluon distribution functions. At the moment we can make an observation that the contribution from minijets to differential elliptic flow is in qualitative agreement with the data .

To understand the centrality dependence of the elliptic flow coming from minijets we use the definition

$$v_2(B)=\frac{d^2p_T𝑑yv_2(p_T,\underset{¯}{B})(dN/d^2p_Tdy)}{d^2p_T𝑑y(dN/d^2p_Tdy)}$$
(43)
to obtain using Eqs. (27) and (31)

$$v_2(B)=\alpha _s\left(\frac{2\mathrm{ln}2\pi K_2N_c}{c^2C_FS_{}(B)Q_s^2}\right)^{1/2},$$
(44)
where $`S_{}(B)`$ is now the $`B`$-dependent nuclear overlap area. Since $`S_{}(B)Q_s^2N_{part}`$ we may therefore conclude that minijets flow scales as

$$v_2(B)\frac{1}{\sqrt{N_{part}}}.$$
(45)
Elliptic flow from Eq. (45) is smaller for central collisions with large $`N_{part}`$ and it increases for peripheral collisions with decreasing $`N_{part}`$, again in qualitative agreement with the data .

## 4 Our model versus experimental data

There are three interesting kinematical regions in high energy heavy ion collisions corresponding to various values of transverse momenta $`p_T`$:

1. $`p_T\genfrac{}{}{0pt}{}{_<}{}Q_s`$ where multiple rescatterings of partons dominate. In this region transverse momentum dependence of $`v_2`$ is given by Eq. (41). Owing to the nonlinear interactions nuclear structure functions in region (i) are functions of only one variable $`Q/Q_s(x)`$ instead of two $`Q`$ and $`x`$. This is often referred to as geometric scaling .
2. $`p_TQ_s`$ where the DGLAP  evolution starts to be important. There the leading behavior of $`v_2`$ is given by Eq. (42) in which DGLAP equation may introduce logarithmic corrections, possibly changing the power of $`\mathrm{ln}(p_T/\mu )`$. The scaling behavior is broken down.
3. At high energies in the course of quantum evolution the saturation scale $`Q_s`$ becomes much larger than any soft QCD scale $`\mathrm{\Lambda }`$. It was argued both analytically  and numerically  that the geometrical scaling holds in a much wider kinematical region than (i). In fact the additional scaling region is $`Q_s\genfrac{}{}{0pt}{}{_<}{}p_T\genfrac{}{}{0pt}{}{_<}{}Q_s^2/\mathrm{\Lambda }`$. Here even though the scattering amplitude is still far from saturation, the nonlinear interactions produce remarkable scaling phenomenon. It was argued in  that the saturation scale in RHIC kinematical region is of the order of $`1÷1.4`$ GeV. Since $`\mathrm{\Lambda }0.1÷0.2`$ GeV, we expect the scaling to be important at $`p_T5÷10`$ GeV, i.e., throughout most of the kinematic region of the differential elliptic flow data .

To compare predictions of our model with data collected by the STAR collaboration at RHIC  we should evaluate $`v_2`$ given by Eq. (40) and Eq. (44) at $`p_T<5`$ GeV. This means calculation in the regions (i) and (iii) where geometric scaling works and the only relevant dimensional scale is $`Q_s`$. Quantum evolution changes the quasi-classical result of Eq. (26) so that the saturation scale acquires energy dependence. Its analytical expression can be calculated from the nonlinear evolution equation . There is a simple way to introduce $`Q_s`$ into our model preserving geometric scaling developed by quantum evolution, where $`Q_s`$ will play a role of a phenomenological parameter which is used to fit the data. Namely, by definition the saturation scale is a scale at which the expression in the Glauber exponent of Eq. (25) equals one

$$\frac{1}{4}\underset{¯}{z}^2Q_s^2(\underset{¯}{z})|_{\underset{¯}{z}^2=4/Q_s^2}=\frac{Q_{s0}^2}{Q_s^2}\mathrm{ln}\frac{Q_s}{2\mu }=\mathrm{\hspace{0.17em}1},$$
(46)
where we have introduced in the quasi-classical (Glauber) approximation

$$Q_{s0}^2=\frac{4\pi \alpha _s^2A}{S_{}^A},$$
(47)
which is an unknown constant for the case of a fully evolved distribution. $`Q_{s0}`$ can be easily expressed in terms of $`Q_s`$ using Eq. (46)

$$Q_{s0}^2=\frac{Q_s^2}{\mathrm{ln}\frac{Q_s}{2\mu }}.$$
(48)
Substituting Eq. (48) into Eq. (40) means the following change in the power of the Glauber exponents

$$\frac{1}{4}\underset{¯}{z}^2Q_{s0}^2\mathrm{ln}\frac{1}{z\mu }=\frac{1}{4}\underset{¯}{z}^2Q_s^2\frac{\mathrm{ln}\frac{1}{z\mu }}{\mathrm{ln}(Q_s/2\mu )}.$$
(49)
Using Eq. (49) in Eq. (40) we can write for the integrals involved

$$_0^{\mathrm{}}\frac{dz}{z^3}J_n(p_Tz)\left(1e^{z^2Q_s^2(z)/4}\right)^2=p_T^2_0^{\mathrm{}}\frac{d\xi }{\xi ^3}J_n(\xi )\left(1e^{(\xi ^2Q_s^2/4p_T^2)\frac{\mathrm{ln}\frac{2p_T}{\xi \mathrm{\Lambda }}}{\mathrm{ln}(Q_s/\mathrm{\Lambda })}}\right)^2$$
(50)
where we defined $`\xi p_Tz`$, $`\mathrm{\Lambda }=2\mu `$ and $`n=0\text{or}\mathrm{\hspace{0.17em}2}`$. Geometric scaling  implies that the distribution functions depend only on one parameter $`p_T/Q_s`$. To eliminate the scaling violating terms in Eq. (50) we note that the average value $`<z>\mathrm{\hspace{0.17em}1}/<p_T>\mathrm{\hspace{0.17em}2}/Q_s`$ and rewrite the logarithm as

$$\mathrm{ln}\frac{2p_T}{\xi \mathrm{\Lambda }}\mathrm{ln}\frac{Q_s}{\xi \mathrm{\Lambda }}$$
(51)
where $`\mathrm{\Lambda }`$ would be an independent parameter of our fit restricted by reasonable possible values of the non-perturbative scale. The integral of Eq. (50) becomes

$$p_T^2_0^{\mathrm{}}\frac{d\xi }{\xi ^3}J_n(\xi )\left(1e^{(\xi ^2Q_s^2/4p_T^2)\frac{\mathrm{ln}\frac{Q_s}{\xi \mathrm{\Lambda }}}{\mathrm{ln}(Q_s/\mathrm{\Lambda })}}\right)^2.$$
(52)
This form of the integrals will be used in Eq. (40) to numerically estimate $`v_2(p_T)`$.

To obtain numerical value of $`v_2(p_T)`$ we have to estimate the normalization coefficient $`K_1`$. To this aim note that the fraction of the total particle multiplicity per unit of rapidity due to soft saturation physics in heavy ion collisions at center-of-mass energy $`\sqrt{s}=130`$ GeV is given by

$$\frac{dN_{sat}}{dy}\mathrm{\hspace{0.17em}1.02}\frac{3N_{\mathrm{part}}}{2},$$
(53)
where $`3/2`$ is the conversion factor from charged particles to all particles multiplicity. On the other hand, we can calculate the total particle multiplicity by integrating Eq. (27) over all momenta $`k_1`$ which gives

$$\frac{dN}{dy}=K_1\frac{2\mathrm{ln}2S_{}^AC_FQ_s^2}{\pi ^2\alpha _s}\mathrm{ln}\frac{Q_s}{2\mu }.$$
(54)
Equating Eq. (53) and Eq. (54) for the head-on collisions ($`B=0`$) we obtain

$$K_1=\frac{3\pi ^2\alpha _s\mathrm{\hspace{0.17em}1.02}N_{\mathrm{part}}(B=0)}{4\mathrm{ln}2S_{}^AC_FQ_s^2\mathrm{ln}(Q_s/\mathrm{\Lambda })}.$$
(55)
where $`N_{\mathrm{part}}(B=0)=344`$ . For gold nuclei with $`\alpha _s=0.3`$, $`Q_s=1.0`$ GeV and $`\mathrm{\Lambda }0.15`$ GeV the normalization factor is $`K_1\mathrm{\hspace{0.17em}0.14}`$.

With the help of Eq. (52) the integration in Eq. (40) has been done numerically for $`\mathrm{\Lambda }0.15`$ GeV. We put $`K_2=2`$ in Eq. (40) in agreement with collinear factorization approaches . $`S_{}^A=\pi R^2`$ with the radius given by Woods–Saxon parameterization $`R=1.1A^{1/3}`$ fm. The results are shown in Fig. 3. It can be seen that

our model describes the high $`p_T`$ data remarkably well. We are led to the conclusion that the high momentum tail of the $`v_2`$ is saturated by the two-particle correlated minijet production. We can also explain the turnover of $`v_2`$ at $`p_T1÷2Q_s`$ GeV as being due to saturation effects in the wave functions of the colliding nuclei at low $`p_T`$. Another important observation is that two-particle correlations in the initial state give significant contribution at low $`p_T`$ as one can see in Fig. 3. Note that the $`p_T^2`$ increase of $`v_2`$ at small momenta may be an artifact of our model which neglects the change of the anomalous dimension $`\gamma `$ of the gluon structure function. It is well known that correct anomalous dimension at the boundary of the saturation region is $`\gamma 1/2`$ . It smoothly changes from 0 in the unitarity limit to 1 in the Bjorken limit . More accurate analysis is needed to understand the low $`p_T`$ behavior of $`v_2`$ due to two-particle correlations in minijet production.

Applicability of our model is restricted to the mid-rapidity kinematical region. However, since the total particle multiplicity is dominated by the particle production at mid-rapidity the averaging over all impact parameters in Eq. (37) gives a reasonable approximation. The situation is different if we want to apply our model to describe the centrality dependence of $`v_2`$. While we expect it to be in agreement with the data for the most central events, it is not legitimate to use it for peripheral collisions.

Eqs. (44) and (45) imply

$$v_2(B)=\frac{v_2(B=0)\sqrt{N_{\mathrm{part}}(B=0)}}{\sqrt{N_{\mathrm{part}}}}=\alpha _s\left(\frac{\pi N_cK_2N_{\mathrm{part}}(B=0)}{2\mathrm{ln}2C_FS_{}^AQ_s^2N_{\mathrm{part}}}\right)^{1/2},$$
(56)
where we used $`c=2\mathrm{ln}2`$ .

In Fig. 4 we show the result of numerical calculations compared to STAR data . The data is represented by black dots while our fit is given by empty crosses. We have used Eq. (56) with the coefficient $`K_2=2`$. To find the average $`N_{\mathrm{part}}`$ in a given centrality bin we have used the model of Kharzeev and Nardi . Horizontal error bars in our fit correspond to the widths of centrality bins. Vertical error bars account for our use of cylindrical nuclear shape instead of spherical and for the uncertainty in $`N_{\mathrm{part}}`$ coming from the model of . As expected our model provides a reasonable description of the data for central events while the fit does not work that well for peripheral collisions. The saturation scale decreases toward the large impact parameters where at some point it becomes of the order of $`\mathrm{\Lambda }_{QCD}`$ and the small coupling approach breaks down. Thus, as was mentioned above, our approach is not applicable to very peripheral collisions and can not be expected to work for them. We conclude that our model gives a reasonable description of $`v_2`$ centrality dependence for the values of centrality $`n_{\mathrm{ch}}/n_{\mathrm{max}}0.4`$.

The maximum of the centrality distribution of elliptic flow is probably determined by the onset of non-perturbative effects. The flow from Eq. (44) increases with increasing $`B`$, until at some point the saturation scale $`Q_s(B)`$ becomes of the order of $`\mathrm{\Lambda }_{QCD}`$ and non-perturbative effects take over keeping $`v_2(B)`$ approximately constant for peripheral collisions. (Strictly speaking Eq. (44) was derived for cylindrical nuclei with saturation scales independent of the impact parameter. However in real life $`Q_s`$ is of course a function of $`B`$ .) Let us denote the impact parameter at which the non-perturbative effects take over by $`B_0`$, so that $`Q_s(B_0)\mathrm{\Lambda }_{QCD}`$. Since $`Q_s`$ is an increasing function of the center of mass energy $`s`$ , the impact parameter $`B_0`$ must also be an increasing function of $`s`$. That is as energy increases more and more collisions become perturbative with the non-perturbative physics being responsible only for increasingly more peripheral collisions. The maximum of the flow centrality distribution could be estimated from Eq. (44) to be $`v_2(B_0)\mathrm{\hspace{0.17em}1}/\sqrt{S_{}(B_0)\mathrm{\Lambda }_{QCD}^2}`$. With the increase of energy $`B_0`$ increases, $`S_{}(B_0)`$ decreases and, therefore, $`v_2(B_0)`$ increases, which qualitatively agrees with RHIC and SPS data .

## 5 Conclusions

In this paper we have calculated the contribution of pairwise azimuthal correlations in minijet production to elliptic flow variable. The resulting differential elliptic flow $`v_2(p_T)`$ given by Eq. (40) is in good qualitative (see Eqs. (41) and (42)) and quantitative (see Fig. 3) agreement with the emerging RHIC data . The centrality dependence of the minijet contribution to elliptic flow $`v_2(B)`$ (see Eqs. (44) and (45)) successfully describes the data for sufficiently central collisions as shown in Fig. 4. The maximum $`v_2(B)`$ appears to increase with energy also in agreement with the data .

There are several important questions which still have to be addressed in the future as more flow data is produced at RHIC. One question concerns the value of the directed flow $`v_1`$. It may seem from the above discussion that using two-particle correlations one might write, similar to Eq. (6),

$$v_1(p_T)=\frac{\mathrm{cos}(\varphi _1(p_T)\varphi _2)}{\sqrt{\mathrm{cos}(\varphi _1\varphi _2)}}$$
(57)
and using the correlations of Eq. (28) get a non-zero directed flow. However, let us remind the reader that in the analysis of directed flow the signs of the weights used in determination of the reaction plane are reversed in backward hemisphere with respect to the forward one . That is, in the numerator of Eq. (57), the contributions from $`y_2>0`$ and $`y_2<0`$ come in with different signs (we denote the central rapidity of the event by $`y=0`$). Therefore, since in our boost-invariant model the two contributions are identical, the final result for $`v_1`$ at mid-rapidity is zero in agreement with the SPS data . If we impose some constraint on the rapidity interval of the analyzed particles making it asymmetric with respect to $`y=0`$, e.g. analyzing only particles with $`y>0`$, the differential directed flow from minijets would also become non-zero as observed experimentally .

In our analysis we have neglected contributions from correlated production of three and more particles in a single subcollision. The corresponding diagrams may also contribute to the two-particle correlation function in Eq. (8). Of course they are suppressed by extra powers of the strong coupling constant $`\alpha _s`$ compared to Eq. (28), though these extra powers are compensated by the large logarithms arising due to phase space integration . In heavy ion collisions at RHIC the saturation scale is of the order of $`Q_s\mathrm{\hspace{0.17em}1}\text{GeV}`$  and the corresponding value of the coupling constant is $`\alpha _s(Q_s)\mathrm{\hspace{0.17em}0.3}`$, which is not too small and will also be enhanced by logarithms of energy and transverse momentum. Thus to have a precise description of the minijet contribution to elliptic flow one has to resum all these multiple emission diagrams, which is equivalent to including the effects of quantum evolution in the double inclusive gluon distribution of Eq. (28). Since at the moment there is even no rigorously derived expression for double gluon production in AA in the quasi-classical approximation, including evolution effects in it seems like an important but not immediate goal.

To improve the quality of the minijet predictions one has to also relax the $`y_1y_2`$ (or $`y_1y_2`$) condition which we imposed to simplify the calculation. This condition led to absence of rapidity correlations in the two-particle distribution given by our Eq. (24). Also the azimuthal correlations given by Eq. (24) are only “back-to-back”, i.e., the produced particles are correlated mostly in the opposite directions in the transverse plane. This behavior seems to contradict recent results reported by PHENIX  where the two-particle correlation function $`C(\mathrm{\Delta }\varphi =\varphi _1\varphi _2)`$ has maxima at both $`\mathrm{\Delta }\varphi =\pi `$ (back-to-back) and $`\mathrm{\Delta }\varphi =0`$. Relaxing the $`y_1y_2`$ (or $`y_1y_2`$) condition would greatly complicate the particle production calculations  though would give predictions for elliptic flow in a much more realistic kinematics. The appropriate calculations were performed for production of a pair of particles with $`y_1y_2`$ by Leonidov and Ostrovsky in . The results of  indicate that the two produced particles are indeed somewhat correlated in rapidity. The analysis of azimuthal distributions performed in  also demonstrates enhancement of both $`\mathrm{\Delta }\varphi =\pi `$ and $`\mathrm{\Delta }\varphi =0`$ correlations in the obtained double particle spectrum, in qualitative agreement with PHENIX data . These correlations should also contribute to smallness of directed flow $`v_1`$. A detailed analysis of elliptic flow employing the full calculation of  will be done elsewhere . In a complete analysis one has to put in realistic geometrical shapes for the nuclei as well. Of course the effects of different cuts used to analyze the data and various choices of weights have to be incorporated in the precise analysis too. While these modifications are likely to introduce some numerical changes in the fits to the elliptic flow data presented above, we believe that the qualitative conclusions of importance of the minijet production to flow variables would remain unchanged.

We therefore conclude that minijets give a large contribution to the elliptic flow extracted using current flow analysis methods, probably accounting for most of the flow at large $`p_T`$ (see Fig. 3). Differential elliptic flow appears to be sensitive to saturation physics of the early stages of the collisions. To study QGP it would be very useful to invent a method of flow analysis that would be insensitive to minijet contribution and would only measure the collective effects due to elliptic flow . We have to note, however, that if one tries to calculate the double inclusive minijet production cross section in the saturation approach  many of the diagrams that would contribute would reduce to independent production of two gluons which would then exchange a gluon with each other. This $`22`$ process would constitute the first step in the onset of thermalization . Therefore it appears that correlated production of pairs of minijets may be intrinsically related to thermalization. On one hand this implies that the minijet flow of Eq. (40) is partially due to collective effects after all. On the other hand it may be impossible to create an observable which would distinguish these early stage thermalization effects from the flow of the fully thermalized quark-gluon plasma.

## Acknowledgements

The authors would like to thank Adrian Dumitru, Miklos Gyulassy, Dima Kharzeev, Roy Lacey, Larry McLerran, Jean-Yves Ollitrault, and Jan Rak for stimulating and informative discussions. Our thanks go to Derek Teaney and Raju Venugopalan for pointing out to us a numerical error in the earlier version of the paper.

The work of Yu. K. was supported in part by the U.S. Department of Energy under Grant No. DE-FG03-97ER41014 and by the BSF grant $`\mathrm{\#}`$ 9800276 with Israeli Science Foundation, founded by the Israeli Academy of Science and Humanities. The work of K. T. was sponsored in part by the U.S. Department of Energy under Grant No. DE-FG03-00ER41132.


text:
# Electromagnetic Transition in Waveguide with Application to Lasers

## Abstract

The electromagnetic transition of two-level atomic systems in a waveguide is calculated. Compared with the result in free space, the spontaneous emission rate decrease because the phase space is smaller, and meanwhile, some resonance appears in some cases. Moreover, the influence of non-uniform electromagnetic field in a waveguide on absorption and stimulated emission is considered. Applying the results to lasers, a method to enhance the laser power is proposed.

PACS codes: 42.50.ct, 42.55.Ah

1 Introduction

The behavior of atoms in confined space had been studied for many years since the possibility of modifying spontaneous emission rates was first mentioned by Purcell . Many studies had been done in theoretical and experimental aspects . The transition rate of atoms may increase  or decrease  in a cave. In this paper, we will discuss spontaneous emission, absorption and stimulated emission of atoms in a waveguide and apply the results to the problem of how to enhance the output power of lasers.

In confined space, the values of momenta are discrete at certain directions and a nonzero lower limit exists due to the boundaries. Compared with free space, the phase space in confined space is smaller. Some processes, e.g. atomic decay, will be influenced. In the case of spontaneous emission, the transition rate of excited atoms is determined by two factors: The transition matrix element between two states and the phase space volume of final state. In confined space, the transition matrix element is the same as in free space under the first order approximation. The phase space volume is smaller because some components of photon momenta in final state are constrained to some discrete values. It leads to the decrease of the total spontaneous emission rate of excited atoms. On the other hand, when the frequency of the photon in final state equals to one of the eigenfrequencies of the cave, resonance will occur and the transition rate will remarkably increase. In this paper, we will give a direct calculation about the spontaneous emission rate of excited atoms in a matrix waveguide. Also, the ratio of the spontaneous emission rate in the waveguide to the case in free space is given. The ratio is determined only by the shape of the boundary, but it is independent of the form of the transition matrix. This means that the ratio is not related to special transition processes. Besides, since the electromagnetic field has a defined non-uniform distribution in a waveguide, the rates of absorption and stimulated emission will depend on the position of the atoms. The corresponding transition rate at each position in the waveguide and the mean transition rate are also given.

The power of lasers is related to the spontaneous-emission lifetime of the excited states. The lifetime of atoms can be influenced in confined space. Therefore, we can change the laser power by putting the atoms in a cavity. In the last part of this paper, we will discuss the change of the output power of lasers under the influence of the scale of waveguides.

2 Spontaneous emission in a waveguide

For one-photon decay of an initial state $`|a`$ into a final state $`|b`$ and a photon with frequency $`\omega _n`$ and polarization $`\sigma `$, the $`S`$-matrix element to first order in perturbation theory is (Here we use Heaviside’s units and take $`c=\mathrm{}=1`$) :

$$S_{ba}=i2\pi \delta (E_b+\omega _nE_a)\frac{1}{\sqrt{L_xL_yL_z}}f_{ba},$$
(1)
where $`E_a`$ and $`E_b`$ are the energies of the two states of the atomic system. $`f_{ba}`$ is the transition amplitude. Under the dipole approximation, the amplitude takes the form:

$$f_{ba}=ie\sqrt{\frac{\omega _n}{2}}ϵ_n^\sigma 𝐫_{ba}.$$
(2)
The vectors $`ϵ`$ are polarization vectors and $`𝐫_{ba}`$ is the matrix element

$$𝐫_{ba}=b|𝐫|a.$$
(3)
Consider a two-level system placed in a waveguide with sides $`L_x`$ and $`L_y`$ which are comparable with the wavelength of the photon. In the waveguide, the components of the momentum at $`x`$ and $`y`$ directions of the photon are constrained to some discrete values:

$$\begin{array}{c}k_x=\frac{n_x\pi }{L_x},k_y=\frac{n_y\pi }{L_y},(n_x,n_y=0,1,2,\mathrm{}).\end{array}$$
(4)
Summing over all final photon states, the total spontaneous emission rate is

$$W_{sp}=\underset{n_xn_yn_z}{}\frac{1}{L_xL_yL_z}2\pi \delta (E_b+\omega _nE_a)|f_{ba}|^2.$$
(5)
Since the photon is free at $`z`$ direction, the summation of $`n_z`$ can be replaced by integration:

$$\begin{array}{ccc}W_{sp}\hfill & =\hfill & \underset{n_xn_y}{}\frac{1}{L_xL_y}\frac{dk_z}{2\pi }2\pi \delta (E_b+\omega _nE_a)|f_{ba}|^2\hfill \\ & =\hfill & \underset{n_xn_y}{}\frac{1}{L_xL_y}\frac{\omega }{k_0}|f_{ba}|^2\hfill \\ & =\hfill & \underset{n_xn_y}{}\frac{1}{L_xL_y}\frac{e^2\omega ^2}{2k_0}|𝐫_{ba}|^2\mathrm{sin}^2\theta ,\hfill \end{array}$$
(6)
where $`k_0`$ is the $`z`$-component of wave vector $`𝐤`$

$$k_0=\sqrt{\omega ^2k_x^2k_y^2}$$
(7)
and $`\theta `$ is the angle between the dipole moment and wave vector $`𝐤`$.

In free space, the total spontaneous emission rate is known as

$$W_{sp}^{free}=\frac{e^2\omega ^3}{3\pi }|𝐫_{ba}|^2.$$
(8)
Comparing with eq.(6), we obtain the ratio of the rates in confined space and free space:

$$\eta =\frac{W_{sp}}{W_{sp}^{free}}=\underset{n_xn_y}{}\frac{1}{L_xL_y}\frac{3\pi }{2\omega k_0}\mathrm{sin}^2\theta .$$
(9)
Because the space in a waveguide is not homogeneous, a new factor $`\mathrm{sin}^2\theta `$ is appeared. This result shows that the spontaneous emission rate is different if the polarization of atoms in the waveguide is different. However, atoms usually are non-polarized, so a mean ratio is needed.

Use ($`\alpha `$, $`\beta `$) to describe the direction of $`𝐫_{ba}`$ in the waveguide. Since

$$\begin{array}{ccc}\mathrm{cos}\theta \hfill & =\hfill & \frac{𝐤𝐞_r}{|𝐤|}\hfill \\ & =\hfill & \frac{k_x\mathrm{sin}\alpha \mathrm{cos}\beta +k_y\mathrm{sin}\alpha \mathrm{sin}\beta +k_z\mathrm{cos}\alpha }{|𝐤|},\hfill \end{array}$$
(10)
the mean ratio can be obtained as

$$\begin{array}{ccc}\overline{\eta }\hfill & =\hfill & \frac{1}{4\pi }\eta 𝑑\mathrm{\Omega }\hfill \\ & =\hfill & \underset{n_xn_y}{}\frac{1}{L_xL_y}\frac{\pi }{\omega k_0}.\hfill \end{array}$$
(11)
Obviously, this result is only related with three length scales: the wavelength of the photon $`\lambda `$, the widths of the waveguide $`L_x`$ and $`L_y`$. Introducing two parameters:

$$\gamma _x=\frac{L_x}{\lambda /2},\gamma _y=\frac{L_y}{\lambda /2},$$
(12)
we can get the following expression:

$$\overline{\eta }=\underset{n_xn_y}{}\frac{1}{\pi }\frac{1}{\sqrt{\gamma _x^2\gamma _y^2\gamma _y^2n_x^2\gamma _x^2n_y^2}}.$$
(13)
This result shows that the spontaneous emission rate is usually smaller in a waveguide than in free space because the phase space of the final state is smaller and the possible states of photon are fewer. However, the spontaneous emission rate will diverge when $`\gamma _x`$ and $`\gamma _y`$ take such values that the denominator of (13) is zero. From (7) we can see that this is just the case that the $`z`$-component of photon $`k_z`$ or $`k_0`$ vanishes. In other words, when $`k_0=0`$, some resonance occurs and the spontaneous emission rate of excited atoms increases rapidly (see fig. 1).

Besides, when both $`\gamma _x`$ and $`\gamma _y`$ are smaller than 1 (see fig. 2), the frequency of the photon is lower than the lowest permitted frequency in such a waveguide. It shows that the phase space of the final state is compressed to zero so the photon can not be emitted and the transition can not occur. At that time, an excited atom will never return to its ground state, i.e., the transition is forbidden.

3 Absorption and stimulated emission in a waveguide

The rates of absorption and stimulated emission in confined space are different from those in free space because the mode structures of electromagnetic fields are different in the two cases. The processes of absorption and stimulated emission are transitions of atomic systems in external electromagnetic fields. The electromagnetic field in a waveguide has a defined distribution which is not uniform, so the transition probabilities will depend on the position of the atoms.

In a matrix waveguide, the expressions for the spatial components of the field are

$$\{\begin{array}{c}E_x(x,y,z)=E_1\mathrm{cos}k_xx\mathrm{sin}k_yye^{ik_zz},\hfill \\ E_y(x,y,z)=E_2\mathrm{sin}k_xx\mathrm{cos}k_yye^{ik_zz},\hfill \\ E_z(x,y,z)=E_3\mathrm{sin}k_xx\mathrm{sin}k_yye^{ik_zz},\hfill \end{array}$$
(14)
where $`k_x,k_y,k_z(=k_0)`$ have been given in (4) and (7), and $`E_1,E_2,E_3`$ satisfy $`k_xE_1+k_yE_2ik_zE_3=0`$.

We will first consider the processes of absorption in the waveguide. In this case, a monochromatic electromagnetic wave at frequency $`\omega `$ is made to interact with an atom which is in state $`|b`$. Since the magnitudes of the atomic electron radii are much smaller than the wavelength of the electromagnetic wave, the dipole approximation is extremely good. On this assumption, we regard $`𝐄`$ as an uniform field in a small region and use its value at $`z=0`$. The corresponding electric field is a sinusoidal function of time $`𝐄(x,y,t)=𝐄_0(x,y)\mathrm{sin}\omega t`$, where $`𝐄_0(x,y)=E_x(x,y,0)𝐢+E_y(x,y,0)𝐣+E_z(x,y,0)𝐤`$. Similar to the case of absorption occurring in free space , the rate of absorption can be obtained directly:

$$W_{ab}(x,y)=\frac{\pi }{3n^2}|𝐫_{ab}|^2\rho (x,y)\delta (\mathrm{\Delta }\omega ).$$
(15)
Here $`n`$ is the refractive index of the atomic system; $`\mathrm{\Delta }\omega `$ is the magnitude of the spread in $`\omega `$; $`𝐫_{ab}`$ is the matrix element corresponding to the transition from the initial state $`|b`$ to the final state $`|a`$; $`\rho `$ is the energy density of the incident electromagnetic wave in the waveguide:

$`\rho (x,y)`$ $`=`$ $`{\displaystyle \frac{1}{2}}n^2E_0(x,y)^2`$
$`=`$ $`{\displaystyle \frac{1}{2}}n^2(E_1^2\mathrm{cos}^2k_xx\mathrm{sin}^2k_yy+E_2^2\mathrm{sin}^2k_xx\mathrm{cos}^2k_yy+E_3^2\mathrm{sin}^2k_xx\mathrm{sin}^2k_yy).`$
Obviously, the transition probability depends on the position of the atom, since the electromagnetic field in waveguide has a non-uniform distribution. The mean rate of absorption in the waveguide is

$$\overline{W}_{ba}=\frac{\pi }{3n^2}|𝐫_{ba}|^2\overline{\rho }\delta (\mathrm{\Delta }\omega ),$$
(17)
with the mean energy density

$$\overline{\rho }=\frac{1}{V}_V\rho (x,y)𝑑V=\frac{n^2}{8}(E_1^2+E_2^2+E_3^2).$$
(18)
The rate of stimulated emission $`W_{ba}`$ is equal to the rate of absorption $`W_{ab}`$, since the transition matrix element $`|𝐫_{ba}|=|𝐫_{ab}|`$ and the photons of stimulated emission must always match the mode of the stimulating photons.

4 Application to lasers

One of the central problems in the area of laser is how to enhance the laser power. There are many factors which can influence the laser power. One of them is the spontaneous-emitted lifetime of excited states. According to the results above, the spontaneous emission rate of excited atoms in confined space, e.g. waveguide, is different to the case of free space. This means that we can change the lifetime of excited atoms by adjusting the shape and scale of the cavity (In this paper, the cavity is a waveguide). One factor which directly influences the laser power is the lifetime of the upper laser level. If we put the laser medium in a waveguide, we can enhance or reduce the laser power by changing the lifetime of the upper laser level.

Take a three-level laser proceeds as an example, the output power of continuous wave laser under the rate-equation approximation is

$$P=A\frac{\chi 1}{\tau }.$$
(19)
Here $`A`$ is a parameter determined by the character of the specific laser medium and laser facility. The quantity $`\tau `$ is the lifetime of the upper laser level, and it is, in general, given by

$$\frac{1}{\tau }=\frac{1}{\tau _{sp}}+\frac{1}{\tau _{nr}},$$
(20)
where $`\tau _{sp}`$ is the spontaneous-emission lifetime and $`\tau _{nr}`$ the nonradiative lifetime of upper laser level. Let $`W_p`$ be the pumping rate and $`W_{cp}`$ the critical pumping rate. The $`\chi `$ is the amount by which threshold is exceeded:

$$\chi =\frac{W_p}{W_{cp}}.$$
(21)
In practice, the following approximation can be taken:

$$W_p\frac{1}{\tau }.$$
(22)
Then the output power is given as

$$P=A(W_p\frac{1}{\tau _{nr}}\frac{1}{\tau _{sp}}).$$
(23)
If the process occurs in a waveguide, $`\tau _{sp}`$ will be replaced by $`\tau _{sp}^W`$ . The $`\tau _{sp}^W`$ is the spontaneous-emission lifetime in the waveguide. From eq.(9), we can obtain

$$\eta =\frac{W_{sp}}{W_{sp}^{free}}=\frac{1/\tau _{sp}^W}{1/\tau _{sp}}=\frac{\tau _{sp}}{\tau _{sp}^W}.$$
(24)
Substitute this relation into eq.(23), the laser power in a waveguide $`P^W`$ is therefore

$$P^W=A(W_p\frac{1}{\tau _{nr}}\frac{\eta }{\tau _{sp}}).$$
(25)
When $`\eta <1`$, the output power enhances, i.e. $`P^W>P`$. On the contrary, when $`\eta >1`$, the output power reduces, i.e. $`P^W<P`$. The analysis carried out here implies that the output power can be adjusted by changing the scale of the waveguide.

Besides, in waveguides, by adjusting the spontaneous-emission lifetime of the upper laser level $`\tau _{sp}^W`$, we can obtain a high power laser pulse. If $`\tau _{sp}^W`$ is long enough, after sufficient pumping, a large population inversion will be obtained. For instance, the spontaneous emission of the upper laser level can be forbidden to a certain extent by modifying the scale of the waveguide. When the population inversion is large enough, change the waveguide to a proper scale to induce a spontaneous transition of the upper laser level. In this way, a laser pulse is obtained.

In conclusion, we have shown that transition rates of atoms may change in confined space. In some cases, the lifetime can be changed by adjusting the experimental arrangements. Based on this idea, a method to enhance the laser power is proposed.

Figure captions:

Fig. 1. In a waveguide, the decay rate is usually depressed since the phase space is smaller than in free space, but remarkably increases when resonance occurs ($`\gamma _y=1.6`$).

Fig. 2. When $`\gamma _x,\gamma _y<1`$, the decay is forbidden (Here $`\gamma _y=0.8`$).