format: text

text:
# Nonlinear waves in a cylindrical Bose-Einstein condensate

## I Introduction

Solitary waves that may occur in a Bose-Einstein Condensate (BEC) have been traditionally discussed in terms of the classical Gross-Pitaevskii (GP) model which is appropriate for the description of weakly correlated systems dalfovo . For instance, a simple soliton was obtained by Tsuzuki tsuzuki  in a homogeneous 1D model, while Zakharov and Shabat zakharov  developed inverse-scattering techniques for the study of multisolitons. Interestingly, the elementary soliton proved to be relevant for an accurate semiclassical description kulish ; ishikawa  of an intriguing mode proposed earlier by Lieb lieb  in a full quantum treatment of a 1D Bose gas based on the Bethe Ansatz lieb2 .

The above developments had long remained purely theoretical because of the absence of a physical realization of a strictly 1D Bose gas. Nevertheless, the picture has significantly changed with the recent observation of similar coherent structures in confined BECs of alkali-metal atoms burger ; denschlag . The very method of experimental production of solitary waves (phase imprinting) was inspired by the analytical structure of the 1D soliton, while various effective 1D models have been developed for their theoretical investigation perez ; jackson ; muryshev ; feder ; dikande ; muryshev2 ; salasnich . On the other hand, the actual stability of the theoretically predicted 1D solitary waves should be questioned within the proper 3D environment of realistic traps muryshev ; feder ; dikande . An important step in that direction was the experimental observation anderson  that a dark soliton initially created in a finite trap eventually decays into vortex rings, as is also predicted by a numerical solution of the corresponding initial-value problem in a 3D classical GP model feder .

Therefore, it is important to carry out a calculation of potential nonlinear modes without a priori assumptions about their effective dimensionality. One could envisage a picture in which the actual solitary waves are hybrids of quasi-1D solitons and 3D vortex rings. It is the aim of the present paper to make the above claim precise by calculating solitary waves that propagate along a cylindrical trap in a steady state with constant velocity $`v`$. Our approach was motivated by the calculation of vortex rings in a homogeneous BEC due to Jones and Roberts jones1  and a similar calculation of semitopological solitons in planar ferromagnets semi .

We have already described the main result of this work in a recent short communication komineas  but a substantial elaboration is necessary in order to appreciate its full significance. Thus the problem is formulated in Sec. II where we also present a brief but complete recalculation of the ground state and the corresponding linear (Bogoliubov) modes for comparison. A detailed calculation of nonlinear modes is given in Sec. III and the main conclusions are summarized in Sec. IV.

## II Formulation and linear modes

The physical picture that we have in mind is a slight idealization of the experiment in Ref. burger . Thus we consider a cigar-shaped trap filled with atoms of mass $`m`$. The transverse confinement frequency is denoted by $`\omega _{}`$ and the corresponding oscillator length by $`a_{}=(\mathrm{}/m\omega _{})^{1/2}`$. The longitudinal confinement frequency $`\omega _{}`$ is assumed to be much smaller than $`\omega _{}`$, hence we make the approximation of an infinitely-long cylindrical trap with $`\omega _{}=0`$. Accordingly, complete specification of the system requires as input the average linear density $`\nu `$ which is the number of atoms per unit length of the cylindrical trap. Finally, we consider the two dimensionless combinations of parameters:

$$\gamma =\nu a,\gamma _{}=\nu a_{},$$
(1)
where $`a`$ is the scattering length related to the coupling constant as usual by $`U_0=4\pi \mathrm{}^2a/m`$.

Now, in the actual experiment of Ref. burger , the trap is filled with 87Rb atoms, the transverse frequency is chosen as $`\omega _{}=2\pi \times 425`$ Hz, the oscillator length is calculated to be $`a_{}0.5\mu m`$, and the estimated linear density lies in the range $`\nu >0.1\text{atoms}/\text{Å}`$. Also taking into account a scattering length $`a50\text{Å}`$, the dimensionless parameters (1) take values in $`\gamma >5`$ and $`\gamma _{}>5\times 10^2`$. In fact, our subsequent calculations will be carried out for a much wider range of the above parameters. Therefore, apart from the idealization of a cylindrical trap, our results are fairly realistic and could be applied to a number of cases of experimental interest.

It is useful to introduce rationalized units through the rescalings

$$t\frac{t}{\omega _{}},\text{r}a_{}\text{r},\mathrm{\Psi }\frac{\nu ^{1/2}}{a_{}}\mathrm{\Psi }.$$
(2)
The energy functional extended to include a chemical potential is then given by

$$W=\frac{1}{2}\left[\mathbf{}\mathrm{\Psi }^{}\mathbf{}\mathrm{\Psi }+\rho ^2\mathrm{\Psi }^{}\mathrm{\Psi }+g(\mathrm{\Psi }^{}\mathrm{\Psi })^22\mu \mathrm{\Psi }^{}\mathrm{\Psi }\right]𝑑V,$$
(3)
where $`g=4\pi \gamma `$ and $`\rho ^2=x^2+y^2`$. Eq. (3) yields energy $`W`$ in units of $`\gamma _{}(\mathrm{}\omega _{})`$ whereas the chemical potential $`\mu `$ is measured in units of $`\mathrm{}\omega _{}`$. The corresponding rationalized equation of motion reads

$$i\frac{\mathrm{\Psi }}{t}=\frac{1}{2}\mathrm{\Delta }\mathrm{\Psi }+\frac{1}{2}\rho ^2\mathrm{\Psi }+g(\mathrm{\Psi }^{}\mathrm{\Psi })\mathrm{\Psi }\mu \mathrm{\Psi },$$
(4)
and depends only on the dimensionless coupling constant $`\gamma `$, because $`g=4\pi \gamma `$ and the chemical potential $`\mu =\mu (\gamma )`$ is fixed by the requirement that the system carry in its ground state a definite average linear density $`\nu `$.

An important first step is thus to obtain accurate information about the ground-state wave function $`\mathrm{\Psi }=\mathrm{\Psi }_0(\rho )`$ which is normalized according to

$$_0^{\mathrm{}}2\pi \rho 𝑑\rho |\mathrm{\Psi }_0|^2=1,$$
(5)
to conform with our choice of rationalized units. The wave function $`\mathrm{\Psi }_0(\rho )`$ is numerically calculated as the minimum of the energy functional, under the constraint (5) that fixes the chemical potential $`\mu (\gamma )`$, by a variant of a relaxation algorithm dalfovo2 . Explicit results are illustrated in Fig. 1 for some typical values of $`\gamma `$ where we also quote the corresponding values of the chemical potential.

The preceding numerical determination of the ground state will provide the basis for all subsequent calculations. However, it is worth mentioning here some limiting cases where the ground state is known analytically. At $`\gamma =0`$,

$$\mathrm{\Psi }_0=\frac{1}{\sqrt{\pi }}e^{\rho ^2/2},$$
(6)
and the chemical potential degenerates to $`\mu =1`$. In the opposite limit, $`\gamma 1`$, one may use the Thomas-Fermi (TF) approximation baym

$$\mathrm{\Psi }_0=\left[\frac{2}{\pi R_{}^2}\left(1\frac{\rho ^2}{R_{}^2}\right)\right]^{1/2}$$
(7)
for $`0\rho R_{}=2\gamma ^{1/4}`$, and $`\mathrm{\Psi }_0=0`$ for $`\rho >R_{}`$. The chemical potential is given accordingly by $`\mu 2\gamma ^{1/2}`$. A comparison with the accurate numerical solution is shown in the inset of Fig. 1 for $`\gamma =100`$. In fact, as we shall see shortly, the TF approximation provides a reasonable description of some quantities of physical interest even for $`\gamma 1`$.

We will also need some information from the linear (Bogoliubov) modes which have already been calculated in the literature to varying degree of completeness zaremba ; stringari ; fedichev . Here we employ a numerical algorithm of our own briefly described as follows. Equation (4) is linearized by inserting $`\mathrm{\Psi }=\mathrm{\Psi }_0+f+i\chi `$ where $`\mathrm{\Psi }_0(\rho )`$ is the calculated ground-state wave function while $`f=f(\text{r},t)`$ and $`\chi =\chi (\text{r},t)`$ are real functions that account for small fluctuations around the ground state. It is somewhat more convenient to use the linear combinations $`a=f+\chi `$ and $`b=f\chi `$ which satisfy the linearized equations

$$\frac{}{t}\left(\begin{array}{c}a\\ b\end{array}\right)=M\left(\begin{array}{c}a\\ b\end{array}\right),M\left(\begin{array}{cc}V& D\\ D& V\end{array}\right)$$
(8)
where

$$V=g\mathrm{\Psi }_0^2,D=\frac{1}{2}\mathrm{\Delta }+\frac{1}{2}\rho ^2+2V\mu .$$
(9)
Our task is then to calculate the spectrum of the differential operator $`M`$ whose eigenvalues are purely imaginary and come in pairs $`\pm i\omega `$ where $`\omega `$ is the sought after physical frequency.

We restrict attention to axially-symmetric waves that propagate along the $`z`$ axis with wave number $`q`$. The Laplace operator is then replaced by

$$\mathrm{\Delta }=\frac{^2}{\rho ^2}+\frac{1}{\rho }\frac{}{\rho }q^2$$
(10)
and the amplitudes $`a`$ and $`b`$ may be assumed to depend only on the radial distance $`\rho `$. A finite-matrix approximation of the operator $`M`$ is obtained by expanding both $`a`$ and $`b`$ in terms of a basis set of non-orthogonal Gaussian wave packets with randomly chosen oscillator lengths papspathis . It is also prudent to enlarge the basis set by including the ground-state wave function $`\mathrm{\Psi }_0(\rho )`$ itself, in order to directly account for the zero (Goldstone) mode associated with the number symmetry. The resulting algorithm is then quite efficient and provides stable approximations of the low-lying eigenvalues even if we include a small number of basis elements.

In Fig. 2 we present explicit results for the lowest eigenfrequency $`\omega =\omega (q)`$ for the same set of coupling constants as in Fig. 1. At $`\gamma =0`$, $`\omega (q)`$ reduces to the free-particle quadratic dispersion $`\omega =q^2/2`$, as expected. At nonzero $`\gamma `$, the dispersion becomes linear near the origin, $`\omega c|q|`$, where $`c`$ is the speed of sound for which explicit values are also quoted in Fig. 2. Finally, we note that our results are in apparent agreement with the Bogoliubov dispersion calculated earlier within the TF approximation zaremba ; stringari  as well as numerically fedichev  – even though a different parameterization of the spectrum was employed in the latter reference.

The speed of sound is a quantity of special physical interest and will also play an important role in the theoretical development of Sec. III. Hence we have carried out a calculation for a wider set of coupling constants and the results are summarized in Fig. 3. It is interesting that our accurate numerical results are consistent with the TF approximation zaremba ; stringari ; kavoulakis

$$c\gamma ^{1/4}$$
(11)
even for values of $`\gamma `$ as low as 1, where the error is about 5%, whereas the error is reduced to less than 1% for $`\gamma \stackrel{>}{}10`$. This fact is especially important because Eq. (11) was employed for the analysis of experimental data andrews . The relative accuracy of this approximation progressively deteriorates in the region $`\gamma <1`$, but a new asymptote, namely

$$c(2\gamma )^{1/2},$$
(12)
was predicted to be reached for sufficiently weak couplings jackson . The weak-coupling approximation (12) is actually consistent with our numerical data for $`\gamma <1/4`$, as is shown in the inset of Fig. 3. However, we should add that the linear part of the Bogoliubov dispersion becomes very narrow in this region of couplings.

## III Nonlinear waves

We now turn to the calculation of axially-symmetric solitary waves traveling along the $`z`$ axis in a steady state with constant velocity $`v`$. These are described by a wave function of the form $`\mathrm{\Psi }=\mathrm{\Psi }(\rho ,\xi )`$, with $`\xi =zvt`$, which is inserted in Eq. (4) to yield the stationary differential equation

$`iv{\displaystyle \frac{\mathrm{\Psi }}{\xi }}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{\Delta }\mathrm{\Psi }+{\displaystyle \frac{1}{2}}\rho ^2\mathrm{\Psi }+g(\mathrm{\Psi }^{}\mathrm{\Psi })\mathrm{\Psi }\mu \mathrm{\Psi },`$ (13)
$`\mathrm{\Delta }`$ $`=`$ $`{\displaystyle \frac{^2}{\rho ^2}}+{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{}{\rho }}+{\displaystyle \frac{^2}{\xi ^2}}.`$
The wave function must vanish in the limit $`\rho \mathrm{}`$, thanks to the transverse confinement, while the condition

$$\underset{\xi \pm \mathrm{}}{lim}|\mathrm{\Psi }(\rho ,\xi )|=|\mathrm{\Psi }_0(\rho )|$$
(14)
enforces the requirement that the local particle density coincide asymptotically with that of the ground state calculated in Sec. II. But the phase of the wave function is not fixed a priori at spatial infinity except for a mild restriction implied by the von Neumann boundary condition

$$\underset{\xi \pm \mathrm{}}{lim}\frac{\mathrm{\Psi }}{\xi }=0$$
(15)
adopted in our numerical calculation. Our task is then to find concrete solutions of Eq. (13) that satisfy the boundary conditions just described.

An important check of the numerical calculation is provided by the virial relation

$$v𝒫=\left[\frac{1}{2}\frac{\mathrm{\Psi }^{}}{\xi }\frac{\mathrm{\Psi }}{\xi }+\rho ^2\mathrm{\Psi }^{}\mathrm{\Psi }+\frac{g}{2}(\mathrm{\Psi }^{}\mathrm{\Psi })^2\mu \mathrm{\Psi }^{}\mathrm{\Psi }\right]𝑑V,$$
(16)
obtained by standard scaling arguments semi . Here $`𝒫`$ is the linear momentum given by the usual definition

$$𝒫=\frac{1}{2i}\left(\mathrm{\Psi }^{}\frac{\mathrm{\Psi }}{z}\frac{\mathrm{\Psi }^{}}{z}\mathrm{\Psi }\right)𝑑V=n\frac{\varphi }{z}𝑑V,$$
(17)
and is measured in units of $`\mathrm{}\nu =\gamma _{}(\mathrm{}/a_{})`$. In the second step of Eq. (17) we employ hydrodynamic variables defined from

$$\mathrm{\Psi }=\sqrt{n}e^{i\varphi },$$
(18)
where $`n=|\mathrm{\Psi }|^2`$ is the local particle density and the phase $`\varphi `$ may be used to construct the velocity field $`\text{u}=\mathbf{}\varphi `$.

Numerical solutions of Eq. (13) are obtained by an iterative Newton-Raphson algorithm jones1 ; semi  briefly described as follows. Suppose that $`\mathrm{\Psi }=\mathrm{\Psi }_{\mathrm{in}}`$ is an initial rough guess for the solution at some velocity $`v`$. We then insert in Eq. (13) the configuration $`\mathrm{\Psi }_{\mathrm{out}}=\mathrm{\Psi }_{\mathrm{in}}+X`$ and keep terms that are at most linear in the amplitude $`X`$. Thus we derive an inhomogeneous differential equation of the form $`LX=Y`$ where the linear operator $`L`$ and the source $`Y`$ are both calculated in terms of $`\mathrm{\Psi }_{\mathrm{in}}`$. We solve this linear system for $`X=L^1Y`$ to obtain $`\mathrm{\Psi }_{\mathrm{out}}=\mathrm{\Psi }_{\mathrm{in}}+L^1Y`$ which is used as input for the next iteration until convergence is achieved at some specified level of accuracy. The procedure is repeated by incrementing the velocity to a different value, typically in steps of $`\delta v=\pm 0.01`$, using as input the converged configuration obtained at the preceding value of the velocity. Therefore, the main numerical burden consists of constructing a finite-matrix lattice approximation of the linear operator $`L`$ which is then inverted by standard routines appropriate for sparse linear systems.

The Newton-Raphson algorithm typically converges after a few iterations and the final configuration is independent of the specific choice $`\mathrm{\Psi }_{\mathrm{in}}`$. But it is also clear that the algorithm will not converge for most choices of $`\mathrm{\Psi }_{\mathrm{in}}`$. Hence it is important to invoke an educated guess for the input configuration provided by the product Ansatz:

$$\mathrm{\Psi }_{\mathrm{in}}=[c_1ic_2\mathrm{tanh}(c_3\xi )]\mathrm{\Psi }_0(\rho )$$
(19)
which capitalizes on the analytically known solitary wave in the homogeneous 1D model tsuzuki ; zakharov  and the ground-state configuration $`\mathrm{\Psi }_0(\rho )`$ numerically calculated in Sec. II. The constants $`c_1,c_2`$ and $`c_3`$ are definite functions of the velocity $`v`$ within the strictly 1D model, but such precise relations need not be invoked for our current purposes except for the normalization condition $`c_1^2+c_2^2=1`$ that is necessary to enforce the boundary condition (14). In other words, the above constants are treated here as trial parameters until we achieve convergence for a specific velocity $`v`$. A corrolary of the preceding discussion is that the converged configuration does not depend on the precise choice of those parameters, and it is certainly not in the form of a product Ansatz often employed for the derivation of effective 1D models jackson ; salasnich . Finally, we note that the Ansatz (19) satisfies the parity relations

$$\text{Re}\mathrm{\Psi }(\rho ,\xi )=\text{Re}\mathrm{\Psi }(\rho ,\xi ),\text{Im}\mathrm{\Psi }(\rho ,\xi )=\text{Im}\mathrm{\Psi }(\rho ,\xi )$$
(20)
which are compatible with Eq. (13) and are actually satisfied by all solutions constructed in the present paper.

We begin with the special case of the relatively weak coupling $`\gamma =1`$ for which the speed of sound was calculated to be $`c=0.95`$ in Sec. II. The simplest possibility is to first attempt to derive a static ($`v=0`$) soliton starting with the input configuration (19) applied for, say, $`c_1=0`$ and $`c_2=1=c_3`$. Indeed, the algorithm quickly converges to a wave function with a nontrivial imaginary part but vanishing real part. The velocity is then incremented to positive values in steps of $`\delta v=0.01`$ and the corresponding wave functions acquire also a nontrivial real part. The process may be continued until the velocity $`v`$ approaches the speed of sound $`c`$ beyond which the solitary wave ceases to exist. An equivalent sequence of solitary waves with velocities in the range $`c<v0`$ is obtained either by starting again with the $`v=0`$ soliton and pushing it to negative velocities or, simply, by taking the complex conjugate of the wave function calculated for $`0v<c`$, since

$$vv,\mathrm{\Psi }\mathrm{\Psi }^{},$$
(21)
is an obvious symmetry of Eq. (13). A detailed illustration of the calculated solitary wave function is given in Fig. 4 for two representative values of the velocity: $`v=c/2`$ and $`v=0`$.

A partial but more transparent illustration is given in Fig. 5 which depicts the level contours of the local particle density $`n=|\mathrm{\Psi }|^2`$ for the two special cases considered in Fig. 4. In words, the calculated solitary wave is a mild soundlike disturbance of the ground state when $`|v|`$ approaches the speed of sound $`c`$, while it becomes an increasingly dark soliton with decreasing $`|v|`$ and reduces to a completely dark (black) soliton at $`v=0`$.

It is now important to calculate the energy-momentum dispersion of the solitary wave. The excitation energy is defined as

$$E=WW_0$$
(22)
where both $`W`$ and $`W_0`$ are calculated from Eq. (3) applied for the solitary wave $`\mathrm{\Psi }(\rho ,\xi )`$ and the ground state $`\mathrm{\Psi }_0(\rho )`$, respectively. The presence of the chemical potential in Eq. (3) provides the compensation that is necessary to compare energies of states with the same number of particles. Similarly, the relevant physical momentum is not the linear momentum $`𝒫`$ of Eq. (17) but the impulse $`Q`$ defined in a manner analogous to the case of a homogeneous gas kulish ; ishikawa ,

$`Q`$ $`=`$ $`{\displaystyle (nn_0)\frac{\varphi }{z}𝑑V}=𝒫\delta \varphi ,`$ (23)
$`\delta \varphi `$ $``$ $`{\displaystyle _0^{\mathrm{}}}2\pi \rho 𝑑\rho n_0(\rho )[\varphi (\rho ,z=\mathrm{})\varphi (\rho ,z=\mathrm{})],`$
where $`n_0=|\mathrm{\Psi }_0(\rho )|^2`$ is the ground-state particle density and $`\delta \varphi `$ is now the weighted average of the phase difference between the two ends of the trap. The delicate distinction between linear momentum and impulse has been the subject of discussion in practically all treatments of classical fluid dynamics batchelor ; saffman  and continues to play an important role in the dynamics of superfluids jones1 . Here we simply postulate the validity of the definition of impulse in Eq. (23) and note that the corresponding group-velocity relation

$$v=\frac{dE}{dQ}$$
(24)
is satisfied to an excellent accuracy in our numerical calculation and thus provides a highly nontrivial check of consistency. In turn, the virial relation (16) is verified using the standard definition of the linear momentum $`𝒫`$ in Eq. (17), as expected. We finally note that the same phase difference $`\delta \varphi `$ which is important for experimental production of solitary waves through phase imprinting burger ; denschlag  is also crucial for the calculation of the impulse.

The dispersion $`E=E(Q)`$ calculated for the complete sequence of solitary waves with velocities in the range $`c<v<c`$ is illustrated in Fig. 6. The apparent $`2\pi `$ periodicity seems surprising, but occured also in the original calculation of a similar mode by Lieb lieb  within a full quantum treatment of a 1D Bose gas interacting via a $`\delta `$-function potential. The Lieb mode was later rederived by a fairly accurate semiclassical approximation based on the elementary solitary wave of the 1D classical GP model kulish ; ishikawa .

Lieb further argued that the corresponding Bogoliubov mode is no more elementary and thus proposed an intriguing dual interpretation of the excitation spectrum. It should be noted that the dispersions of the two modes exhibit the same linear dependence at low momenta, $`Ec|Q|`$, where $`c`$ is the speed of sound, but significant differences arise at finite momenta. The differences are especially pronounced in the current calculation within a cylindrical trap. Specifically, let us assume an average linear density $`\nu =0.02\text{atoms}/\text{Å}`$ which leads to $`\gamma =\nu a=1`$ and $`\gamma _{}=\nu a_{}=10^2`$. If we then adjust the $`\gamma =1`$ Bogoliubov dispersion of Fig. 2 to the units employed in Fig. 6, the two dispersions are seen to diverge very quickly at the scale of Fig. 6. In other words, Bogoliubov and Lieb modes operate at rather different energy and momentum scales in a realistic trap.

To summarize the preceding accurate calculation for $`\gamma =1`$, the solitary wave is essentially quasi-1D in this weak-coupling region and its main features are indeed captured by an effective 1D model jackson2 . However, all cases of actual experimental interest burger ; denschlag ; anderson  are characterized by significantly larger values of the effective coupling where quasi-1D solitons are expected to be unstable muryshev ; feder . In particular, Ref. muryshev  suggests that a critical coupling occurs in a cylindrical trap when $`n_{\mathrm{max}}U_0/\mathrm{}\omega _{}=2.4`$, where $`n_{\mathrm{max}}`$ is the maximum local particle density in the ground state of the trap. If we tentatively assume that the TF approximation (7) can be trusted at the anticipated critical coupling, the above criticality condition reads $`2\gamma ^{1/2}=2.4`$ where $`\gamma `$ is the dimensionless effective coupling constant defined in Eq. (1). The root of this equation $`\gamma =\gamma _c=1.44`$ provides a critical coupling $`\gamma _c`$ above which quasi-1D solitons are predicted to be unstable jackson2 .

In fact, our numerical calculation suggests that the quasi-1D nature of the solitary wave is completely lost at a higher critical coupling, namely for $`\gamma >\gamma _13.9`$. The emerging new picture is clear at $`\gamma =10`$ which is the special case described in our recent short communication komineas . This case is reanalyzed and further extended in the continuation of the present paper.

It is natural to begin again with the calculation of a static ($`v=0`$) soliton obtained by using the input configuration (19) with $`c_1=0,c_2=1`$ and practically any $`c_3`$. We then increment the velocity to both positive and negative values in steps of $`\delta v=\pm 0.01`$ to yield a sequence of solitary waves which now display two surprising features. First, a ringlike structure develops for $`\gamma =10`$ that was not present at $`\gamma =1`$. Second, the above sequence exists only over the limited velocity range $`v_1<v<v_1`$ where $`v_1=0.84=0.47c`$ and $`c=1.77`$ is the speed of sound calculated in Sec. II for $`\gamma =10`$. The existence of a critical velocity $`v_1`$ also becomes apparent in the energy-momentum dispersion of the above sequence depicted by a dotted line in Fig. 7. This portion of the dispersion is symmetric around $`Q=\pi `$, where it achieves a maximum, but remains open ended at two critical points that correspond to $`v=\pm v_1`$.

It is thus not surprising that an independent sequence of solitary waves with lower energy exists for $`\gamma =10`$. Indeed, we return to the input configuration (19) but now target a solution with velocity in the range $`v_1<v<c`$. After some experimentation a solution is obtained for, say, $`v=1.5`$ if we choose the trial parameters $`c_1=0.2`$, $`c_2=0.98`$ and $`c_3=3`$. Having thus obtained a specific solution for $`v=1.5`$ the algorithm is iterated forward and backwards in steps of $`\delta v=\pm 0.01`$ to obtain an entirely new sequence of solitary waves in the velocity range $`v_1<v<c`$, and a corresponding sequence for $`c<v<v_1`$ through the symmetry relation (21). Here $`v_1`$ is the same critical velocity encountered in the preceding paragraph, as is also apparent in the calculated energy-momentum dispersions which are depicted by solid lines in Fig. 7 and join the previously calculated dotted line through cusps that correspond to $`v=\pm v_1`$.

Hence we turn to a description of the detailed nature of this new sequence of solitary waves. For values of the velocity near the speed of sound $`c`$, the calculated soliton appears again as a mild soundlike disturbance of the ground state. The dominant features of the solitary wave are pronounced as the velocity is decreased to lower values and become reasonably apparent for $`v=c/2`$ that corresponds to point I in the dispersion of Fig.7. The wave function is completely illustrated through its real and imaginary parts in Fig. 8. An important new feature emerges by comparison with the corresponding case at $`\gamma =1`$ illustrated in frame I of Fig. 4. Both the real and the imaginary part at the center of the soliton ($`\xi =0`$) now vanish for a specific radius $`R=2.8`$, thus a vortex ring is beginning to emerge. A partial but more transparent illustration is given in Fig. 9 where we depict the radial dependence of the local particle density $`n=|\mathrm{\Psi }|^2`$ for various values of $`\xi `$. Again it is clear that the density near the center of the soliton ($`\xi =0`$) vanishes on a ring with a relatively large radius $`R=2.8`$. The features of the vortex ring become completely apparent, and its radius is tightened, as we proceed to smaller values of the velocity. A notable special case is the static ($`v=0`$) vortex ring with radius $`R=1.8`$ illustrated in frame II of Fig. 9, which is far from being a black soliton. The corresponding point II in Fig. 7 is thus a new local maximum of the energy-momentum dispersion, which is clearly distinguished from the local maximum at point IV that corresponds to the static black soliton discussed earlier in the text.

One would think that pushing the velocity $`v`$ to negative values would somehow retrace the calculated sequence of vortex rings backwards. In fact, our algorithm continues to converge to vortex rings of smaller radii until the critical velocity $`v_1`$ is encountered where the ring achieves its minimum radius $`R_{\mathrm{min}}=0.8`$ and ceases to exist for smaller values of $`v`$. The terminal state at $`v=v_1`$ is illustrated in frame III of Fig. 9. We have thus described a sequence of solitary waves that consists of bonafide 3D vortex rings and does not contain a black soliton. The corresponding branch in the energy-momentum dispersion of Fig. 7 is labeled by points I, II, and III that stand for the special cases $`v=c/2,0`$, and $`v_1`$. As mentioned already, an equivalent sequence of solitary waves exists in the range $`c<v<v_1`$ and leads to a dispersion curve in Fig. 7 that is mirror symmetric to the branch (I,II,III) around $`Q=\pi `$.

To complete the description for $`\gamma =10`$ we must briefly return to the auxiliary sequence of solitary waves associated with the portion of the dispersion that is depicted by a dotted line in Fig. 7. As one moves from point III to point IV, the ringlike structure is more or less preserved at constant radius $`R=R_{\mathrm{min}}=0.8`$. Nevertheless, the detailed features of the vortex ring are tamed at small velocities and completely disappear for $`v=0`$ to yield a black soliton at point IV.

We thus essentially conclude our description of solitary waves for $`\gamma =10`$ by schematically summarizing our main results in Fig. 10. Yet some of the elements of the preceding discussion are sufficiently surprising to deserve closer attention. For example, simple inspection of Fig. 7 reveals that the group velocity becomes negative in the region (II,III) or, equivalently, the impulse is opposite to the group velocity. This rotonlike behavior is consistent with the Onsager-Feynman view of a roton as the ghost of a vanished vortex ring donnelly  because the calculated radius of the vortex ring is monotonically decreasing along the fundamental (I,II,III) sequence. A full-scale roton would develop if the terminal point III were an inflection point beyond which the group velocity begins to rise again. Actually, this is exactly what happens as one moves away from point III along the upper branch in Fig. 7, but this “roton” portion of the dispersion now appears in a strange location by comparison to the usual situation in liquid helium donnelly . On the other hand, the black soliton at the stationary point IV is indeed the ghost of a vanished vortex ring, as explained in the preceding paragraph.

It is also interesting to question how the picture just described evolves with increasing values of the dimensionless coupling constant $`\gamma `$ which is the only parameter that enters the rationalized GP equation. Our numerical calculations have revealed yet another critical coupling $`\gamma _212`$, in the sense that new flavor arises for $`\gamma >\gamma _2`$. The structure of the solitary waves in this new regime becomes sufficiently clear for $`\gamma =20`$ and is best summarized by the calculated energy-momentum dispersion shown in Fig. 11. Apart from mirror symmetry, the dispersion now exhibits two cusps that correspond to two critical velocities $`v_1=1.35=0.64c`$ and $`v_2=0.48=0.23c`$, where $`c=2.1`$ is the speed of sound calculated for $`\gamma =20`$ as in Sec. II.

The nature of the solitary waves associated with the various branches in the dispersion of Fig. 11 is very briefly described with the aid of Fig. 12. Thus we consider the sequence of five characteristic points (I,II,…,V) that roughly cover half of the dispersion, the other half being obtained by the mirror symmetry (21). The lowest branch (I,II) corresponds to single vortex rings with velocities in the range $`c>v>v_1`$, as in the case $`\gamma =10`$. Again the ring achieves its minimum radius at the critical velocity $`v=v_1`$ (point II). The new element for $`\gamma =20`$ is the intermediate branch (II,III,IV) that corresponds to double rings with velocities in the range $`v_2>v>v_1`$. The second ring is first created at the flanks of the trap and comes closest to the original ring at the new critical velocity $`v=v_2`$ (point IV). This double-ring configuration is more or less preserved along the upper branch (IV,V), with velocities in the range $`v_2>v>0`$, but gradually fades away to become a black soliton at $`v=0`$ (point V).

While the numerical calculation becomes increasingly more difficult for larger values of $`\gamma `$, it is clear that a sequence of critical couplings $`\gamma _1,\gamma _2,\mathrm{}`$ exists and leads to a hierarchy of axisymmetric vortex rings. The single-ring solution associated with the lowest branch of the spectrum is a robust feature for all $`\gamma >\gamma _13.9`$ and is likely stable to all perturbations. But it is possible that the multiple-ring configurations associated with the higher branches are unstable to non-axisymmetric perturbations. In this respect, one should recall that the solitary wave that corresponds to the upper branch in the original calculation of Jones and Roberts jones1  within the homogeneous GP model was later argued to be unstable jones2 .

However, we should emphasize that the vortex rings constructed here differ significantly from the Jones-Roberts (JR) vortex ring that provided the basic motivation for the present work. As with ordinary smoke rings in fluid dynamics, the JR ring can never be static thanks to a virial relation of the type (16) that prevents finite-energy solutions with $`v=0`$ in the homogeneous GP model. As a result, the radius of the vortex ring grows to infinity at low velocity. This picture is completely rearranged in a cylindrical trap because the occurrence of slow vortex rings with large radius is restricted by the boundaries of the trap. Instead, vortex rings are predicted to nucleate at the flanks of the trap as soundlike pulses with high velocity approaching the speed of sound $`c`$, and their radius actually decreases with decreasing velocity. In particular, it is now possible to obtain static ($`v=0`$) vortex rings of finite radius that are no longer contradicted by the virial relation (16). The structure of the energy-momentum dispersions calculated throughout the present paper clearly reflects the substantial restructuring of vortex rings within a cylindrical trap.

It is then natural to question whether or not there exists a limit in which the JR vortex ring is recovered. One should expect that this may happen when the bulk healing length $`L=(na)^{1/2}`$ is significantly smaller than the transverse oscillator length $`a_{}`$. This limit is translated into large values of $`\gamma =\nu a`$ which is the only dimensionless parameter that enters the rationalized GP equation. Now, our discussion earlier in this section suggests an increasingly complicated hierarchical structure in the strong-coupling limit rather that a simple JR soliton. A logical conclusion is that a JR vortex ring somehow created within the bulk will sooner or later sense the boundaries of the cylindrical trap. It will thus either directly dissipate into sound waves, or reorganize itself to conform with one or more of the presently calculated vortex rings possibly after ejecting some amount of radiation in the form of sound waves.

The preceding remarks indicate a certain non-uniformity that is inherent in the approximation of the cigar-shaped trap by an infinitely-long cylindrical trap. The same phenomenon is also apparent in the calculation of linear modes in Sec. II. For instance, neither one of the two asymptotes for the speed of sound quoted in Eqs. (11) and (12) approaches the well-known Bogoliubov speed in a homogeneous Bose gas jackson . But a non-uniformity of this type is not a reason to doubt that a sufficiently elongated trap can be approximated by an infinitely-long cylindrical trap.

## IV Conclusion

We have thus presented a complete 3D calculation of solitary waves in a cylindrical Bose-Einstein condensate. In all cases considered there exists a nontrivial phase difference $`\delta \varphi `$ that is reminiscent of strictly 1D solitons tsuzuki ; zakharov  and is important for their experimental production through phase imprinting burger ; denschlag .

Nevertheless, the detailed structure of the solitary waves depends crucially on the strength of the dimensionless effective coupling constant $`\gamma `$. Quasi-1D solitons occur only in the weak-coupling region $`\gamma <\gamma _13.9`$ where some of our accurate numerical results could be approximated through effective 1D models jackson2 . But a sufficiently strong coupling or density is necessary in order to pronounce the special features of a condensate. It is thus not surprising that the effective coupling in experiments performed so far lies in the region $`\gamma >\gamma _1`$ where the nature of the theoretically predicted solitary waves changes drastically.

For $`\gamma >\gamma _1`$ solitary waves are still characterized by a nontrivial phase difference $`\delta \varphi `$ between the two ends of the trap but are otherwise 3D vortex rings. This finding is consistent with a recent experiment anderson  and a corresponding theoretical analysis feder  in finite traps. The main mathematical advantage of the approximation of a sufficiently elongated trap by an infinitely-long cylindrical trap is that vortex rings can then be calculated in a steady state propagating with a constant velocity $`v`$. It is thus possible to carry out a detailed study of the soliton profile as a function of the effective coupling constant $`\gamma `$ and the velocity $`v`$, as is done in the present paper.

An interesting by-product of the above idealization is that a soliton is characterized by a definite energy-momentum dispersion. The calculated dispersion is found to be the direct analog of the Lieb mode lieb  in the weak-coupling region and acquires interesting rotonlike features for stronger couplings. Perhaps such a dispersion can be measured by a combination of phase imprinting burger ; denschlag  with Bragg spectroscopy recently employed for the detection of the usual Bogoliubov mode stamper ; ozeri . In this respect, one should keep in mind that Bogoliubov and Lieb modes operate at different energy and momentum scales in a realistic trap. This fact becomes evident by the different sets of physical units employed for the Bogoliubov mode in Fig. 2 and the Lieb mode in, say, Fig. 6. The difference is accounted for by the second dimensionless coupling $`\gamma _{}=\nu a_{}`$ in Eq. (1) which is much stronger than $`\gamma =\nu a`$ because $`\gamma _{}/\gamma =a_{}/a10^2`$.

###### Acknowledgements.

We are grateful to A.R. Bishop, N.R. Cooper, G.M. Kavoulakis, F.G. Mertens, and X. Zotos for valuable comments.


text:
# References

NOTES ON RENORMALIZATION

A. A. Vladimirov

BLTP, JINR, Dubna, Moscow region 141980, Russia
alvladim@thsun1.jinr.ru
http://thsun1.jinr.ru/~alvladim/qft.html

We outline the proofs of several principal statements in conventional renormalization theory. This may be of some use in the light of new trends and new techniques (Hopf algebras, etc.) recently introduced in the field.

The following statements form a backbone of renormalization theory  :

* In renormalizable models, all ultraviolet divergences are removed by the use of a special recursive subtraction procedure ($`R`$-operation) .
* This operation is completely equivalent to adding the properly chosen local singular counterterms to the Lagrangian.
* A freedom in the renormalization procedure is precisely a possibility to change parameters in the Lagrangian, via the renormalization group transformations.
* This freedom (renormalization-scheme dependence) cancels out entirely in observable quantities, if expressed in proper variables.

Our aim here is to recall generic proofs for these statements. This may prove to be useful, especially when comparing conventional methods with the new techniques (Hopf algebras, etc.) recently proposed . Addressing the conceptual matter only, we choose to deal with the simplest cases where, nevertheless, the principal ideas of renormalization can be demonstrated in unabridged form. For example, one can imagine a massless theory with a single coupling constant $`g`$ , like $`g\phi ^3`$ in 6 dimensions, additionally ‘quenched’ by discarding all self-energy (propagator-type) subgraphs. Of course, such a model can be non-unitary as quantum field theory, but it is fairly well renormalizable, with only 3-point graphs being divergent, and without overlapping divergences.

We begin with the statement ($`𝐑`$) . Recursive nature of the $`R`$-operation is well understood from the following picture (empty circle designates the *full* 3-*vertex*, i.e., formal sum of all graphs having three external legs, with proper combinatoric factors) :

$$R\text{}=(1K)\text{},\text{}R\text{},$$
(1)
with $`K`$ being a projector, $`K^2=K`$ , picking up the singular (and $`1K`$ the regular) terms in the corresponding regularization parameter. In a word, to renormalize a graph one has to renormalize its subgraphs and then subtract the remaining superficial singularity. Of course, this procedure terminates for any given graph. But there arises a nontrivial question: how to make sure that a singular contribution we thus obtain is local, i.e., does not depend on (or is at most polynomial in) external momenta? This is crucial for renormalizability, because subtracting nonlocal quantities cannot be interpreted (and so justified) as adding counterterms to Lagrangian. As a rule, the proof of locality heavily depends on the regularization prescription used, and, thus, not so much can be said in a general language we are using here. However, the following trick often works: try to find an operation $`D`$ (say, differentiation with respect to external momenta ) which eliminates superficial divergences and commutes with $`K`$ . Then $`KD`$ yields zero when applied to a graph which diverges only superficially, whereas $`DK`$ differentiates the corresponding counterterm, with zero result, thus demonstrating its locality.

Let us now proceed with (C) . The basic reason why recursive subtractions are equivalent to counterterms lies in the auto-recursive nature of the perturbation series itself: an insertion of the full 3-vertex instead of some elementary vertex in a given 3-point graph produces (a subset of) the graphs of the same full 3-vertex, with correct combinatoric coefficients. Namely, the following identity holds:

$$\text{}=\text{}+\text{}+\frac{1}{2}\text{}+\mathrm{}$$
(2)
Here each (except first) term of the sum in the r.h.s. is a ‘primitive’ (only superficially divergent) graph with all its vertices replaced by full 3-vertices. The sets of graphs present in both sides of (2) evidently coincide. The combinatoric factors also agree, due to the same symmetry features of the elementary and the full 3-vertex.

We now *define* a (minimally) renormalized 3-vertex as follows: the coupling constant $`g`$ , at each its occurence in the perturbative expansion of (unrenormalized) full 3-vertex, should be replaced by a series $`g_B=g+𝒪(g^3)`$ (denoted by a black circle below) with the terms $`𝒪(g^3)`$ assumed to be singular in a regularization parameter:

$$\text{}\text{}\text{}\text{}R\text{}$$
(3)
By our assumption,

$$K\text{}=\text{}\text{},K\text{}=0.$$
(4)
Thus, we lay counterterms in the basis of renormalization. Now our goal is twofold: to present an explicit formula for counterterms, and to relate the recipe (3) with $`R`$-operation. Both aims are achieved simultaneously: eq. (2), being identical in $`g`$ , remains an identity after the formal substitution $`gg_B`$ , so

$$\text{}=\text{}+\text{}+\frac{1}{2}\text{}+\mathrm{}$$
(5)
Using (4) immediately yields

$$\begin{array}{c}\text{}=\text{}K\left(\text{}+\frac{1}{2}\text{}+\mathrm{}\right)\hfill \\ \hfill \text{}+S_K\left(\text{}+\frac{1}{2}\text{}+\mathrm{}\right),\end{array}$$
(6)
where the notation $`S_K`$ is used for counterterms. At last, from (6) and (5) we conclude:

$$\begin{array}{c}\text{}=\text{}+(1K)\left(\text{}+\frac{1}{2}\text{}+\mathrm{}\right)\hfill \\ \hfill R\text{}=\text{}+R\left(\text{}+\frac{1}{2}\text{}+\mathrm{}\right),\end{array}$$
(7)
$$S_K\left(\text{}+\mathrm{}\right)=K\left(\text{}+\mathrm{}\right).$$
(8)
Now (6), (8) give a recipe for evaluating counterterms and (7) for renormalizing arbitrary graphs (or, rather, the interplay of these relations provides an iterative procedure for finding both counterterms and renormalized values). Of course, these are readily recognized as the well-known recursive prescriptions of $`R`$-operation in its conventional form.

Now we are in a position to justify the statement (G) . Since (the properly organized) subtractions are done by counterterms, any possible finite (re)subtractions can also be done by finite (non-singular) counterterms. Thus, for observable quantities, all the difference between two renormalization schemes boils down to a (finite) redefinition of $`g`$ ,

$$\stackrel{~}{f}(\frac{p}{\mu },\stackrel{~}{g}(g))=f(\frac{p}{\mu },g),\stackrel{~}{g}(g)=g+𝒪(g^3),$$
(9)
where $`p`$ represents momenta, and $`\mu `$ is a *renormalization parameter* inevitably present in any renormalization procedure (and denoted here by the same letter in both schemes). Another form of this relation addresses the case of different $`\mu `$’s within the same scheme (for different $`\mu `$’s imply different renormalization) :

$$f(\frac{p}{\mu },g)=f(\frac{p}{\mu _0},g_0),g=\overline{g}(\frac{\mu }{\mu _0},g_0).$$
(10)
This property is called *renormalization invariance* . The corresponding transformations $`\mu _0\mu ,g_0g`$ , which leave physics intact, form the *renormalization group* . A function $`\overline{g}`$ is *effective* (or *running*) coupling constant. It shows how $`g`$ must change to compensate a change in $`\mu `$ , and enters into the main formula of the renormalization group approach,

$$f(\frac{p}{\mu },g)=f(1,\overline{g}(\frac{p}{\mu },g)),$$
(11)
which trades the $`p`$-dependence of the function $`f`$ for its dependence on the (now effectively $`p`$-dependent) coupling $`\overline{g}`$ . Furthermore, being actually a record of the group transformation law, $`\overline{g}`$ is subject to standard conditions

$$\overline{g}(xy,g)=\overline{g}(x,\overline{g}(y,g)),\overline{g}(1,g)=g,$$
(12)
which are usually cast into differential equations governed by $`\beta `$-function:

$$(x_x\beta (g)_g)\overline{g}(x,g)=0,x_x\overline{g}(x,g)=\beta (\overline{g}(x,g)),\beta (g)_x\overline{g}(x,g)_{x=1}$$
(13)
The very existence of a (non-singular) dependence $`g(\mu )`$ given by (10) which leaves observable quantities invariant, has interesting consequences . Since both $`g`$ and $`\mu `$ come into play through $`g_B`$ only, one can find the derivative (at $`g_B=\text{const}`$) $`dg/d\mathrm{log}\mu =\beta (g)`$ directly from the (singular) expansion of $`g_B`$ in powers of a regularization parameter $`M`$ :

$$g_B(\frac{M}{\mu },g)=g+a(g)\mathrm{log}\frac{M}{\mu }+b(g)\mathrm{log}^2\frac{M}{\mu }+\mathrm{},$$
(14)
$$\frac{dg}{d\mathrm{log}\mu }=\frac{g_B/\mathrm{log}\mu }{g_B/g}=\frac{a(g)+2b(g)\mathrm{log}\frac{M}{\mu }+\mathrm{}}{1+a^{}(g)\mathrm{log}\frac{M}{\mu }+\mathrm{}}=a(g)+(2b(g)a^{}(g))\mathrm{log}\frac{M}{\mu }+\mathrm{}$$

The first term in the r.h.s. gives $`\beta (g)=a(g)`$ , while other terms are singular and must cancel out. Thus we obtain $`2b(g)=a^{}(g)`$ and other formulas expressing the coefficient functions at higher powers of $`\mathrm{log}(M/\mu )`$ in (14) in terms of $`a(g)`$ .

Finally, consider the statement (S) . From (9) and (11) one can see that the renormalization scheme dependence of $`f`$ and $`\overline{g}`$ is governed by the same function $`\stackrel{~}{g}(g)`$. To show that the two dependences actually cancel, we define (up to a constant) a function $`\psi (g)`$ ,

$$\psi ^{}(g)\frac{1}{\beta (g)}d\mathrm{log}\mu =d\psi (g)\mathrm{log}\mu =\psi (g)+\text{const}$$
(15)
that allows us to rewrite the original definition of $`\overline{g}`$ in the form

$$\psi (\overline{g}(\frac{\mu }{\mu _0},g_0))\psi (g_0)=\mathrm{log}\frac{\mu }{\mu _0}$$
(16)
and introduce (again up to a constant) a new expansion parameter $`\mathrm{log}\mathrm{\Lambda }`$ :

$$\overline{g}(\frac{p}{\mu },g)=\psi ^1(\mathrm{log}\frac{p}{\mu }+\psi (g))=\psi ^1(\mathrm{log}\frac{p}{\mathrm{\Lambda }}),\mathrm{\Lambda }\mu e^{\psi (g)}.$$
(17)
The scheme dependence of $`\psi `$ , according to the rule adopted in (9) : $`\mu \mu ,g\stackrel{~}{g}(g)`$ , is extracted from the last equality in (15) :

$$\stackrel{~}{\psi }(\stackrel{~}{g}(g))=\psi (g)+\text{const},\stackrel{~}{\psi }^1(x+\text{const})=\stackrel{~}{g}(\psi ^1(x)).$$
(18)
Now we come to a scheme-independent expansion recipe :

$$f(\frac{p}{\mu },g)=f(1,\overline{g}(\frac{p}{\mu },g))=f(1,\psi ^1(\mathrm{log}\frac{p}{\mathrm{\Lambda }}))=\mathrm{\Phi }(\mathrm{log}\frac{p}{\mathrm{\Lambda }}),$$
(19)
$$\mathrm{\Phi }(x)f(1,\psi ^1(x)),\stackrel{~}{\mathrm{\Phi }}(x+\text{const})=\mathrm{\Phi }(x).$$
(20)
Let us summarize the whole procedure. One obtains $`f`$ and $`\beta (g)`$ by evaluating Feynman graphs; $`\psi (g)`$ is found from (15) . The scheme dependences of $`f`$ (9) and $`\psi `$ (18) cancel each other in the function $`\mathrm{\Phi }`$ , which is expanded in $`\mathrm{log}(p/\mathrm{\Lambda })`$ like this:

$$\mathrm{\Phi }(\mathrm{log}\frac{p}{\mathrm{\Lambda }})=1+\frac{a_1}{\mathrm{log}\frac{p}{\mathrm{\Lambda }}}+\frac{a_2+a_3\mathrm{log}(\mathrm{log}\frac{p}{\mathrm{\Lambda }})}{\mathrm{log}^2\frac{p}{\mathrm{\Lambda }}}+𝒪(\frac{1}{\mathrm{log}^3\frac{p}{\mathrm{\Lambda }}}).$$
(21)
Numerical coefficients $`a_i`$ may depend on the scheme used, but in a trivial way: all the difference between $`\{a_i\}`$ and $`\{\stackrel{~}{a}_i\}`$ is eliminated by a proper rescaling $`\mathrm{\Lambda }c\mathrm{\Lambda }`$ , with $`c`$ constant. For instance, imposing $`a_2=0`$ renders all other coefficients scheme-independent. In (21), $`\mathrm{\Lambda }`$ is treated as a formal expansion variable, the same for each scheme. However, (17) and (18) explain why $`\mathrm{\Lambda }`$ is suited here: by definition, it behaves like an (almost) scheme-independent quantity, $`\stackrel{~}{\mathrm{\Lambda }}=\mu \mathrm{exp}(\stackrel{~}{\psi }(\stackrel{~}{g}(g)))=\mu \mathrm{exp}(\psi (g)+\text{const})=c\mathrm{\Lambda }`$ .

Acknowledgements. This work was supported by the RFBR grant 00-01-00299.


text:
# Transient X-ray Binaries in Elliptical Galaxies

## 1. Introduction

Chandra  observations of nearby elliptical and S0 galaxies (e.g. Sarazin, Irwin & Bregman 2000; Blanton, Sarazin & Irwin 2000; Sarazin, Irwin, & Bregman 2001; Finoguenov & Jones 2001; Angelini et al. 2001; Kraft et al. 2001) often find $``$ 20 X-ray sources with luminosities greater than $`10^{38}`$ erg sec-1and even more with luminosities greater than $`10^{37}`$ erg sec-1. For a neutron star (NS) accretor, this corresponds to accretion rates as large as $`10^8M_{}`$ yr-1, implying a short lifetime compared to the age of the stellar population. One way to generate such rapid accretion rates in an old stellar population (neglecting formation of tight binaries in dense environments, such as globular clusters) is to have a Roche lobe-filling binary driven by nuclear expansion of a low-mass star on the red giant branch (orbital periods are days or longer). Ritter (1999) has shown an analytic solution for such a Low Mass X-Ray Binary (LMXB), which is further simplified in ellipticals where the initial companion masses should all be $`0.9M_{}`$ for an old population.

However, such wide binaries are subject to thermal instabilities in the accretion disk (King, Kolb & Buderi 1996; see Dubus, Hameury & Lasota 2001 for an up-to-date discussion) that cause transient accretion. Using the current instability criterion, we follow the evolution of such an LMXB and find that over 75% of its life is spent as a transient. We thus expect that a large fraction of the field LMXBs in elliptical galaxies are long-term transients. Our ability to detect them with repeated Chandra observations depends on their outburst durations, which can be as long as years.

We overview the binary evolution scenario in §2 and summarize the challenge of finding the bright, blue optical counterparts to these LMXBs in elliptical galaxies. In §3, we apply the current understanding of accretion disk instabilities in X-ray sources and estimate the outburst recurrence times. In §4, we compare our work to galactic systems, noting that the recently discovered galactic transients with outburst durations in excess of ten years (e.g. Wijnands 2002) might well be wide binaries. We close in §5 with our conclusions and a discussion of alternative binary scenarios that lead to more stable mass transfer, especially relevant for those LMXBs found in globular clusters (Sarazin, Irwin & Bregman 2001; Angelini, Loewenstein & Mushotzky 2001).

## 2. Binary Evolution with Red Giant Donors

We only consider binaries with initial orbital periods $``$1 day. If the binary is tighter ($`P_{\mathrm{orb}}<1`$ day), rapid accretion may be initiated due to orbital angular momentum loss from magnetic braking, resulting in mass transfer while a low-mass star is still on the main sequence (see Tauris & Savonije 1999). Whenever there are multiple parameters to choose, we always take the value that maximizes the stability of the accretion disk. This assures that our calculated fraction of the evolutionary lifetime spent as a transient is a firm lower limit.

We use Ritter’s (1999) analytical solution for the LMXB evolution, which assumes that the donor star is Roche lobe filling on the red giant branch and expanding due to an increasing luminosity. The mass transfer is slow enough that the donor star always remains close to thermal equilibrium; thus $`R_2=R_{2,R}`$, where $`R_2`$ is the radius of the donor star and $`R_{2,R}`$ is the radius of its Roche lobe. For simplicity, we assume that mass and angular momentum are conserved, which implies that at the end of the LMXB’s evolution the NS (of mass $`M_1`$) will have a mass of 1.7-2.1$`M_{}`$. We assume that the initial donor mass is $`0.9M_{}`$ with a population I metallicity, consistent with recent determinations of the age and metallicity of the predominant stellar populations in these elliptical galaxies (Trager et. al. 2000).

Once the evolution is solved, all parameters of the binary can be found in terms of the orbital period at the onset of mass transfer, $`P_i`$, which is measured in days, and the donor’s mass, $`M_2`$, which is measured in solar masses. The mass transfer rate within the binary is

$`\dot{M_2}=6.4\times 10^{10}M_{}\text{ yr}^1\text{ }P_i^{14/15}\left({\displaystyle \frac{M_2}{0.9}}\right)^{7/3}`$
$`\times \left({\displaystyle \frac{2.3M_2}{1.4}}\right)^{14/5}\left({\displaystyle \frac{5}{3}}{\displaystyle \frac{1}{M_2}}{\displaystyle \frac{2}{2.3M_2}}\right)^1,`$ (1)
where $`1.4M_{}`$ is the initial NS mass and $`2.3M_{}`$ is the total mass of the binary. Since we will show in §3 that the disk is likely unstable, we will not concern ourselves here with the implied super Eddington accretion rates at large $`P_i`$ should this mass transfer rate proceed onto the NS. Using equation (1) the evolution of the binary system is shown for varying $`P_i`$ by the solid lines in Figure 1. The long dashed line is the boundary where the hydrogen envelope of the donor has vanished, halting the evolution at the ending orbital period noted. The evolving orbital period is

$`P_{\mathrm{orb}}=P_i\left({\displaystyle \frac{M_2}{0.9}}\right)^3\left({\displaystyle \frac{2.3M_2}{1.4}}\right)^3,`$ (2)
shown in Figure 2 by the solid lines.

For such an LMXB in an elliptical galaxy, we expect: (1) the orbit must be wide in order to reach such high $`\dot{M}`$’s, and (2) the reprocessing of the X-rays by the large accretion disk should make it a bright, very blue object. At the Eddington luminosity, an X-ray binary with $`P_{\mathrm{orb}}=10`$ days (much like Cygnus X-2) has $`M_V=2`$ (van Paradijs & McClintock 1994) and is blue ($`BV0,UB1`$; van Paradijs & McClintock 1995). For NGC 4697 (at a distance modulus of $`30.35`$; Tonry et al. 2001), this corresponds to an apparent V magnitude of 28.4.

Assuming that the LMXB is not associated with a globular cluster, the brightest star which could confuse a search for the optical counterpart would be a red giant branch star just about to reach helium core ignition (a TRGB star). For metallicities $`[\mathrm{Fe}/\mathrm{H}]<1`$, $`M_V(TRGB)2.5`$, whereas for $`[\mathrm{Fe}/\mathrm{H}]0.7`$, $`M_V(TRGB)1.5`$ (Saviane et al. 2000). For a slightly metal-rich star (such as is likely the case for these ellipticals; Trager et al. 2000), the peak brightness is even fainter, $`M_V(TRGB)1.0`$ (Garnavich et al. 1994). Hence, the optical counterpart to a wide LMXB will be brighter at $`V`$ (and especially so in $`B`$ or $`U`$) than any red giant branch star in an elliptical galaxy.

The surface brightness of an elliptical galaxy would make such a search difficult in the central regions. However, beyond $``$ 2 arcminutes, the surface brightness at $`V`$ (in NGC 4697, for example) is below the sky level (for HST) of 22-23 $`\mathrm{Mag}/\mathrm{arc}\mathrm{sec}^2`$ (Goudfrooij et al. 1984). In that limit, and using the measured Chandra position of sources at such locations, such a wide (10 day) binary would be detectable at B or V in a $`5\times 10^4\mathrm{s}`$ observation with the Advanced Camera for Surveys on HST. Of course, such optical counterparts are easy to find in the nearby bulge of M31, where $`V22`$.

## 3. Disk Instability, Transient Fraction, and Recurrence Times

If the disk is thermally stable and transfers matter to the NS at the supplied rate, there is nearly a one-to-one relation between $`P_i`$ and the X-ray luminosity (see Figure 1) that would allow the orbital period distribution to be inferred from the luminosity distribution (Webbink et al. 1983). However, this argument fails because these wide binaries are almost always transient accretors.

The usual condition used for the disk instability is that some portion of the disk is below the hydrogen ionization temperature ($`6500`$ K). Since the effective temperature of the disk decreases with increasing radius, it is usually adequate to check that the outer radius of the disk is below the critical temperature. In X-ray binaries, however, this temperature is fixed by irradiation from the central accreting source (van Paradijs 1996). The surface temperature is then estimated as (Dubus et al. 2001)

$`T_{\mathrm{irr}}^4{\displaystyle \frac{C\eta \dot{M_1}c^2}{4\pi \sigma R^2}},`$ (3)
where $`R`$ is the disk radius, $`\sigma `$ is the Stefan-Boltzmann constant, $`\eta =0.1`$ is the conversion efficiency for accretion onto the NS, and $`C`$ measures the fraction of the X-ray luminosity that reaches the disk (and contains information on the irradiation geometry). We take $`C=5.7\times 10^3`$ (slightly larger than the preferred value, $`C=5\times 10^3`$ of Dubus et al. 2001) so as to maximize the persistent fraction. The disk is assumed to extend out to 70% of the Roche lobe around the NS (King et al. 1997).

Setting $`T_{\mathrm{irr}}=6500`$ K allows us to solve for the critical time, accretion rate, and orbital period at which the instability sets in as a function of the donor’s mass. The critical mass transfer rate is then

$`\dot{M}_{2,\mathrm{crit}}=3.3\times 10^9M_{}\text{ yr}^1M_2^{14/9}\left(2.3M_2\right)^{14/9}`$
$`\times \left({\displaystyle \frac{5}{3}}{\displaystyle \frac{1}{M_2}}{\displaystyle \frac{2}{2.3M_2}}\right)^{10/3},`$ (4)
which is used for the instability boundary shown in Figures 1 and 2 as the dotted line. The critical orbital period is

$`P_{\mathrm{orb},\mathrm{crit}}=5.5\text{ days }M_2^{7/6}\left(2.3M_2\right)^{5/3}`$
$`\times \left({\displaystyle \frac{5}{3}}{\displaystyle \frac{1}{M_2}}{\displaystyle \frac{2}{2.3M_2}}\right)^{5/2},`$ (5)
and this is used for the instability boundary in Figure 2. The fraction of time spent transient is then found (numerically to better than a percent) as a function of the initial orbital period,

$`\mathrm{Transient}\mathrm{Fraction}=0.0500\chi ^2+0.2385\chi +0.7655,`$ (6)
where $`\chi =\mathrm{log}_{10}\left(P_i/\mathrm{d}\right)`$. For $`P_i`$ in excess of a day, the binary is transient at least $`75\%`$ of the time, and this does not include binaries containing black holes which are transient almost all the time (e.g. King et al. 1996).

We estimate a maximum recurrence time for outbursts by assuming the outside-in model. We assume that the binary is quiescent until a critical surface density in the outer accretion disk (at radius $`R`$) has been reached (Dubus, Hameury,& Lasota 2001),

$`\mathrm{\Sigma }_{\mathrm{max}}=644\text{ g cm}^2(M_1/M_{})^{0.37}(R/R_{})^{1.11},`$ (7)
where we assume that the viscosity parameter $`\alpha =0.1`$ and that there is no irradiation of the disk in quiescence. This expression gives a maximum recurrence time

$`t_{\mathrm{rec},\mathrm{max}}=15\text{ yrs }P_i^{1.14}M_2^{3.88}\left(2.3M_2\right)^{2.75}`$
$`\times \left({\displaystyle \frac{5}{3}}{\displaystyle \frac{1}{M_2}}{\displaystyle \frac{2}{2.3M_2}}\right),`$ (8)
that should be viewed with some caution given the uncertainties in the accretion disk physics. This recurrence time is $``$ 100-10,000 years (depending on the initial orbital period; see Figure 2) and is consistent with current estimates from observations of galactic ms radio pulsar/He white dwarf binaries (Ritter & King 2002). Our estimate is maximal because it is possible (and maybe even likely) that the LMXB may outburst before the outer disk reaches the critical mass, such as an outburst initiated in the inner disk.

## 4. Comparison with Galactic Sources

One difficulty with comparing these evolutionary models with observations is the lack of long orbital period LMXBs in our galaxy. Many persistent sources have low enough orbital periods ($`120`$ hours) that they are probably not accurately described by the evolutionary model used above. However, a lack of persistent sources at $`P_{\mathrm{orb}}>1`$ day is in agreement with our calculations and generally expected from previous work (for example King, Frank, Kolb & Ritter 1997). Figure 1 shows that such NS binaries spend most of their lifetime as transient rather than persistent sources ($`>75\%`$ of the time). Though all these relations are consistent, it is still hard to claim that the instability criterion has been fully scrutinized by the galactic population.

Indeed, the outburst durations of these systems might be so long that they appear to us as “persistent” sources (with accretion rates most likely in excess of the mass transfer rate expected from binary evolution) until they suddenly disappear. Indeed, one such source (and maybe others; Wijnands 2002) has now been clearly identified: KS 1731-260 (Wijnands et al. 2001; Rutledge et al. 2002). It was accreting at $`\dot{M}3\times 10^9M_{}\mathrm{yr}^1`$ for at least a decade prior to its sudden decline in luminosity by a factor of $`>10^4`$ (Rutledge et al. 2002). Analysis of the immediate quiescent emission in terms of thermal emission from the NS surface yielded constraints on the recurrence time from 300-1,000 years (Wijnands et al. 2001; Rutledge et al. 2002). Until this system’s orbital period is found (current knowledge of the infrared counterpart points to a possibly evolved companion; Orosz, Bailyn & Whitman 2001) we will not be able to say whether it is an LMXB of the type we are discussing. However, its existence does point to a likelihood of finding such similar systems in elliptical galaxies.

One galactic binary that does have a long orbital period is GRO J1744-28 (11.8 days; Finger et al. 1996). This system is transient and was shown by Bildsten and Brown (1996) to have a companion with mass $`M_20.334M_{}`$ and core mass $`M_c0.22M_{}`$ (if the companion is a population I metallicity star). When this binary’s evolutionary path is predicted and plotted, it is found to be transient (as it is) with a maximum outburst recurrence time of $``$1,000 yrs (see dot-dashed line in Figure 2), consistent with a single outburst.

## 5. Conclusions and Discussion

We have shown that the expected mass transfer scenario for low mass X-ray binaries in the field of an elliptical galaxy (where the donors must be $`<M_{}`$) that yield $`\dot{M}10^8M_{}\mathrm{yr}^1`$ (as observed by Chandra) have thermally unstable disks. Such binaries are expected to be transient accretors with maximum recurrence times of 100-10,000 years for at least $`75\%`$ of their lifetimes. When active, it appears likely that they would easily reach the NS’s Eddington limit.

Our inability to calculate the duration of the bright outbursting phase makes it difficult to predict how many such transients will be identified in repeated Chandra visits. If the duty cycle is bright for ten years out of a thousand (like implied for KS 1731-260; Rutledge et al. 2002), and there are 100 such binaries, then a year later, roughly 10 would have faded into quiescence to be replaced by new transients. This is not that far off from what was seen over a five month interval in Cen A (Kraft et al. 2001), where four sources (out of 120 exposed in both epochs) were either non-detectable ($`L_x<(12)\times 10^{36}\mathrm{erg}\mathrm{s}^1`$) or appeared bright ($`L_x>4\times 10^{37}\mathrm{erg}\mathrm{s}^1`$) in either epoch. There were a larger number ($`35`$) which showed significant variability (Kraft et al. 2001).

Such comparisons are also complicated by the discovery that many of the X-ray sources in ellipticals are coincident with GCs (Sarazin et al. 2001; Angelini et al. 2001; White, Sarazin & Kulkarni 2002). Other possible accretion scenarios then include accretion from lower-mass main sequence stars that are formed in this dense stellar environment (see Rasio, Pfahl & Rappaport 2000 for a recent discussion). These tight systems would most likely be persistent accretors, thus apparently reducing the fraction of transients amongst the detected sources. This might just be an impression however, as if the duty cycle of the field transients is really 1%, then the underlying population is indeed larger than that in GCs. Namely, it is possible that the GC population is only apparently  a large fraction because LMXBs in GCs are persistently bright. The census of the much more numerous population of field transients can only be taken over a few thousand years.

We thank G. Dubus, J. Irwin, H. Ritter, and C. Sarazin for helpful discussions and the anonymous referee for a very informative critique of our work. This research was supported by NASA via grant NAG 5-8658 and by the NSF under Grants PHY99-07949 and AST01-96422. L. B. is a Cottrell Scholar of the Research Corporation.


text:
# Super-poissonian photon statistics and correlations between pump and probe fields in Electromagnetically Induced Transparency

## Abstract

We have measured the photon statistics of pump and probe beams after interaction with Rb atoms in a situation of Electromagnetically Induced Transparency. Both fields present super-poissonian statistics and their intensities become correlated, in good qualitative agreement with theoretical predictions in which both fields are treated quantum-mechanically. The intensity correlations measured are a first step towards the observation of entanglement between the fields.

Electromagnetically Induced Transparency (EIT) is an interference effect that can be observed when three-level atoms interact simultaneously with two lasers. The recent interest in these systems has been driven by observations of very slow light pulse propagation in EIT media  and even of light storage . Even though it is normally viewed as a quantum interference effect, EIT has a very simple classical counterpart . The question of whether there exist intrinsic quantum effects (with no classical analog) in EIT naturally arises in this context. It is even more important in connection with suggested applications in the field of quantum information . Quantum field properties in EIT have been theoretically investigated. Phase-noise squeezing was predicted by Fleischhauer and co-workers . Agarwal  found matched photon-statistics for two classical and two quantum fields interacting with three-level atoms, in a situation of Coherent Population Trapping. Jain  extended this work, and predicted excess-noise correlations in EIT. In all the above , at least one of the fields was treated classically. As we will see below, it is our belief that such an approximation is not valid in the EIT situation.

In this Letter we present the first experimental investigation (to our knowledge) of field fluctuations and correlations in EIT. We have performed photon statistics measurements of both pump and probe fields, as a function of the probe field detuning. As predicted by our theoretical treatment, both fields, treated quantum-mechanically, have super-poissonian statistics. Furthermore, they are coupled by their interaction with the atoms and thus their intensities become correlated.

We begin by describing our theoretical model for three-level atoms in the $`\mathrm{\Lambda }`$-configuration (ground states $`|1`$ and $`|2`$, and excited state $`|0`$) interacting with two quantum modes of the electromagnetic field. In order to simulate the interaction of the atoms with two propagating fields, we use the input-output formalism , and consider the interaction between the two fields and the atoms in a ring cavity , with external input fields. For the intracavity field operators we used plane wave and quasi-monochromatic approximations. The interaction hamiltonian is obtained with the usual dipole and rotating-wave approximations

$$\widehat{H}_{\text{int}}=\mathrm{}g_1\widehat{S}_1^+(t)\widehat{A}_1(t)+\mathrm{}g_2\widehat{S}_2^+(t)\widehat{A}_2(t)+\text{h.c.},$$
(1)
where $`\widehat{A}_1`$ ($`\widehat{A}_2`$) is the annihilation operator for intracavity field 1 - pump (intracavity field 2 - probe), $`g_1`$ ($`g_2`$) is the atom – field 1 (field 2) coupling strength, and $`\widehat{S}_1^+`$ ($`\widehat{S}_2^+`$) the atomic polarization on the transition $`|1|0`$ ($`|2|0`$).

From the Heisenberg equations of motion, we derive quantum Langevin equations  for the system operators. In a matrix form, we find

$$\frac{d\widehat{\overline{X}}(t)}{dt}=\widehat{\text{A}}(t)\widehat{\overline{X}}(t)+\text{D}\widehat{\overline{}}(t).$$
(2)
The vector operator $`\widehat{\overline{X}}(t)`$ has as its components the operators giving the atomic inversions and polarizations (and their hermitian conjugates) corresponding to the two transitions, the coherences between the two ground states and the annihilation and creation operators of both fields. The vector of the Langevin forces acting on the system is given by $`\widehat{\overline{}}(t)`$, and D is the diffusion matrix. The matrix $`\widehat{\text{A}}(t)`$, in the steady state, will represent the drift matrix. Eq. (2) is a compact form of writing twelve coupled differential equations for the operators.

We are interested in studying the fluctuations of the field operators. This is done by linearizing the operators in Eq. (2) around their stationary values. A new set of twelve coupled differential equations is obtained. It is interesting to examine the equations for the fluctuations of one field and for the corresponding atomic polarization:

$$\frac{d\delta \widehat{A}_1}{dt}=\left(\frac{\gamma }{2}+i\mathrm{\Delta }_{c1}\right)\delta \widehat{A}_1i\frac{g_1}{\tau }\delta \widehat{S}_1^{}+\sqrt{\frac{\gamma }{\tau }}\delta \widehat{A}_{1in},$$
(3)
$`{\displaystyle \frac{d\delta \widehat{S}_1^{}}{dt}}`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{\Gamma }_1}{4}}i\delta _{L1}\right)\delta \widehat{S}_1^{}+ig_1w_1\delta \widehat{A}_1+ig_1\alpha _1\delta \widehat{W}_1`$ (5)
$`ig_2s_{12}^{}\delta \widehat{A}_2ig_2\alpha _2\delta \widehat{S}_{12}^++\widehat{F}_{S1}.`$
Here we define $`\gamma `$ as the cavity linewidth, $`\mathrm{\Delta }_{c1}`$ cavity detuning for field 1, $`\tau `$ cavity length divided by the speed of light, $`\widehat{A}_{1in}`$ the annihilation operator for the input field 1, $`\mathrm{\Gamma }_1`$ spontaneous emission rate from $`|0|1`$, $`\delta _{L1}`$ detuning between field 1 and the corresponding atomic transition, $`w_1`$ steady-state inversion between states $`|0`$ and $`|1`$, $`\alpha _1`$ ($`\alpha _2`$) steady-state amplitude of field 1 (field 2), $`\widehat{W}_1`$ inversion (operator) between states $`|0`$ and $`|1`$, $`s_{12}^{}`$ steady-state coherence between ground states $`|1`$ and $`|2`$, $`\widehat{S}_{12}^+`$ coherence operator, $`\widehat{F}_{S1}`$ Langevin fluctuation force. The notation $`\delta \widehat{S}_1^{}`$ means fluctuations of the corresponding operator.

It is clear from Eqs. (3) and (5) that the fluctuations in field 1 are not only determined by input field fluctuations but also by atomic fluctuations and by fluctuations in field 2. Moreover, in Eq. (5) we notice that noise correlations between the fields will arise, conditioned to the existence of a coherence between the ground states, $`s_{12}^{}`$, only important on EIT resonance.

The terms $`ig_1w_1\delta \widehat{A}_1`$, $`ig_1\alpha _1\delta \widehat{W}_1`$, $`ig_2s_{12}^{}\delta \widehat{A}_2`$, $`ig_2\alpha _2\delta \widehat{S}_{12}^+`$ in Eq. (5) are all of the same order. This means that we can not neglect the fluctuations of either field with respect to its average (steady-state value). It is not valid to treat one field as classical and the other as quantized. This conclusion also applies to the work by Agarwal  and Jain  who consider four independent fields (two classical and two quantized) interacting with coherent population trapped (CPT) atoms. The fields used to produce CPT also have inherent quantum fluctuations, which lead to atomic fluctuations and, therefore, should not be neglected.

From the twelve coupled differential equations for the fluctuations, we calculate the spectral matrix S$`(\mathrm{\Omega })`$  for the noise, which yields the fluctuations of the output cavity fields, after Fourier transform of the input-output relations.

The use of a cavity to simulate an EIT experiment with cw propagating fields (no cavity) is adequate as long as we restrict ourselves to fluctuations around the steady state and do not investigate transient effects. We must also consider the cavity large enough so that the atoms experience neither significant changes in their spontaneous emission rates nor collective effects.

The noise spectra calculated for probe and pump amplitude quadratures and their correlations are plotted in Fig. 1, as a function of the probe frequency. These spectra were calculated for an analysis frequency $`\mathrm{\Omega }=\mathrm{\Gamma }/6\pi `$ (consistent with the experiment), where $`\mathrm{\Gamma }`$ is the total spontaneous emission rate from the excited state. We took equal coupling constants for both transitions, pump strictly on resonance, and a ratio of intensities \[$`|\alpha _1/\alpha _2|^2=9`$ in Fig. 1 (a)\] such that there is a very narrow and deep EIT resonance. For the correlations in Fig. 1 (b), both fields have equal intensities. The number of atoms was of the order of $`10^4`$, which is sufficiently large to justify the linearization method used.

We find super-poissonian statistics for both fields, and intensity correlations peaked at the EIT resonance. Physically we understand this behavior as follows. The coupling of the fields with the atoms introduces a spreading in the photon distributions associated to each field. The atoms act as “beamsplitters”, redistributing photons between both modes, in such a way that the mean numbers do not vary appreciably but the variances are modified. Taking absorption into account, we expect this effect to be reduced, since absorption introduces randomness in the photon redistribution process. We therefore expect to observe a maximum effect for minimum absorption, which corresponds to the EIT resonance. The predicted correlations also depend on the presence of absorption. In the EIT situation, the transmission of one field depends on the (fluctuating) intensity of the other field (cross Kerr effect), leading to intensity correlations. The experimental results presented below are in very good qualitative agreement with these predictions.

Our experimental setup is sketched in Fig. 2. Measurements were performed in a 5 cm long room-temperature Rb vapor cell. Two independent extended-cavity diode lasers (ECDLs) provide the pump and probe fields. A small portion of each beam is extracted and sent to auxiliary saturated absorption cells, used as frequency references. The pump beam is tuned to the 85Rb $`|5S_{1/2},F=3|5P_{3/2},F^{}=3`$ transition, while the probe beam is scanned across the Doppler-broadened $`|5S_{1/2},F=2|5P_{3/2},F^{}=1,2,\text{or }3`$ line. The two beams have orthogonal polarizations and are combined by means of a polarizing beamsplitter cube (PBS). At the output of the cell we can separate the beams again with another PBS. We then have the option of detecting either beam with a balanced detection setup . Half-wave plates also enable us to recombine the two fields and measure the photon statistics of the sum and difference intensities of the two beams. From these we can extract the intensity correlations.

The usual EIT signals are observed by sending only the probe beam into the detection region and measuring the DC (average) intensity. The intensity transmitted through the cell has a narrow peak as a function of the probe frequency. The width is much narrower than the Doppler width, which is a signature of the interference effect. Next, we measure the photon statistics (intensity fluctuations) for the probe beam, yielding the signal presented in Fig. 3 (a). For the range of frequencies spanned, the photon statistics is always super-poissonian, but it presents a sharp peak (almost 20 dB) corresponding to the EIT resonance. In this measurement, the initial photon statistics of the probe beam is poissonian (for the frequency of analysis chosen, the pump initially has super-poissonian statistics, but the theory predicts this effect even when both fields have poissonian statistics). The pump photon statistics presents similar behavior, as seen in Fig. 3 (b). In Fig. 3 (c) and (d) we present the corresponding Fano factors for both probe and pump fields, respectively. The Fano factor is given by the ratio of the intensity fluctuations to the average intensity (shot noise level), and we plot it on a linear scale. In these measurements, the probe and pump intensities were 14.5 mW/cm2 and 63.8 mW/cm2, respectively. The observed behavior agrees very well with the theoretical predictions of Fig. 1.

The noise spectra were recorded using a spectrum analyzer (HP 8560 E) in the zero span mode ($`RBW=300`$ kHz, $`VBW=3`$ kHz), with a center (analysis) frequency of 2 MHz. This frequency is in a window such that the probe laser has poissonian photon statistics and the electronic noise is safely below shot noise (more than 5 dB). The pump laser’s noise is $`8`$ dB above shot noise (and shows little variation with the analysis frequency). The analyzer is triggered by the same signal used to scan the probe frequency. The photodetectors are EG&G FND-100 with a nominal efficiency of $`70\%`$. We checked that none of the measurements is influenced by saturation in the detection. This is done by introducing neutral density filters before detection and observing that the photon statistics becomes poissonian, with a linear dependence on absorption.

For typical EIT signals, in the presence of a strong pump field, the shot noise level also peaks, as expected, following the increasing DC intensity on EIT resonance. However, if we lower the pump power, this effect tends to disappear, since incoherent effects then dominate over the coherent pumping on the transition $`|5S_{1/2},F=3|5P_{3/2},F^{}=3`$. The probe DC signal, and consequently the shot noise level, no longer show any evidence of a coherent effect in the atomic medium for pump intensities lower than 0.45 mW/cm2, for a probe intensity of 13.3 mW/cm2 (it is the ratio of pump to probe intensities that is relevant). On the other hand, by looking at the photon statistics, the coherent effect is still clearly identifiable.

We have also measured the noise properties of the sum and difference beam intensities. This is done by sending each beam to one of the detectors. By subtracting the intensity-difference noise from the intensity-sum, we observe a clear correlation on the EIT resonance. The corresponding shot-noise level is measured by mixing the two beams so that each detector receives half of each beam. The results are presented in Fig. 4. We observe $`7`$ dB splitting between the sum and difference fluctuations. For photon statistics, this corresponds to a correlation term $`2C=4\{n_1n_2n_1n_2\}`$.

In summary, we have measured the photon statistics of probe and pump beams after interaction with an atomic medium, in a situation of Electromagnetically Induced Transparency. We observe super-poissonian statistics, peaked on EIT resonance. We have also observed intensity correlations between the two initially independent fields, as a result of their interaction with the atoms. The experimental result features agree very well with our theoretical predictions. The balanced detection scheme leads to enhanced sensitivity for detecting coherent effects in the atomic medium. The importance of the correlations predicted and observed is two-fold. On the one hand, this is a very strong indication that, when investigating quantum properties of one of the fields in a coherent situation (such as EIT), one can not treat either field as classical, as done by several authors. The coherent exchange of photons between the fields, which is at the very origin of the EIT effect, introduces correlations between their fluctuations. Another important issue is the nature of the correlations between both fields. Quantum correlations (entanglement) between intense fields can be used in the context of quantum information. This would be the dual of the experiment by Julsgaard and co-workers , in which correlations between atoms in two different vapor cells were created using light. In the present experiment, we have made the first step in this direction by demonstrating the existence of correlations. However, in order to distinguish between classical and quantum correlations, we will need an improved experimental setup. A new calculation shows that entanglement between pump and probe fields can be produced in the EIT situation, and its observation will require measurements of quadrature fluctuations of both fields .

We thank Profs. L. M. Narducci and H. M. Nussenzveig for critically reading the manuscript, and S. Simionatto for technical assistance. One of us (M.F.S.) would like to thank Prof. S. Salinas for his hospitality in São Paulo. This research was funded by FAPESP. Additional support was provided by the Brazilian agencies CAPES and CNPq.