Resonant processes and exciton propagation in a frozen gas

J. S. Frasier^{∗} and V. Celli^{†}

Department of Physics, University of Virginia, Charlottesville, Virginia 22904 (September 12, 2000)

Abstract

We present numerical simulations of a theory of resonant processes in a frozen gas of excited atoms interacting via dipole-dipole potentials that vary as r^{−3}, where r is the interatomic separation. The simulations calculate timedependent averages of transition amplitudes and transition probabilities for a single atom in a given state interacting resonantly with a uniformly distributed random gas of atoms in a different state. The averages are over spatial configurations of the gas atoms, which are held fixed while the resonant interaction creates a Frenkel exciton that can travel from atom to atom. We check that the simulations reproduce previously known exact results when the exciton is not allowed to propagate [Phys. Rev. A 59, 4358 (1999)]. Further, we develop an approximation for the average transition amplitude that compares well with the numerical results for a wide range of values of the system parameters.

PACS numbers: 34.10.+x, 34.90.+q

Typeset using REVTEX

∗Electronic address: jsf9k@virginia.edu

†Electronic address: vc@virginia.edu

I. INTRODUCTION

Frozen gases are a new and in some ways ideal laboratory to test and refine our understanding of quantum theory in a complex system. With laser technology, the translational temperature of the gas can be lowered to the point where it can be ignored when discussing electronic processes, and at the same time one can manipulate and detect these processes in the gas with an extraordinary selectivity and precision. Frozen Rydberg gases have the added advantage that their states are well understood and electronic processes occur on a microsecond time scale, easily allowing for time-resolved spectroscopy. These gases provide an excellent means of studying many-body effects, because the dipole-dipole interactions which are important for binary collisions at room temperature are still dominant in the frozen gas, but now involve many atoms that are, to a first approximation, standing still.

Pioneering experiments on resonant processes in frozen Rydberg gases have been carried out by Anderson et al. [1,2] and Mourachko et al. [3]. The work presented here was motivated, in part, by the desire to understand those experiments and the subsequent work of Lowell et al. [4] highlighting the dynamic aspects of these resonant, many particle systems. The measurements of Refs. [1,2,4] are of mixtures of ^{85}Rb atoms initially prepared in the 23s and 33s states, henceforth to be called the s and s′ states, respectively. There are initially N atoms in the s state and N′ atoms in the s′ state. The transition ss′ → pp′ is monitored, where p refers to a 24p_{1/2} state and p′ refers to a 34p_{3/2} state, at and near resonance. What is measured is N_{p′}/N′, where N_{p′} is the number of atoms in the p′ state after an elapsed time t. (Measuring N_{p}/N gives nothing really new, since N_{p} = N_{p′} at all times.) Some of the experimental features we wish to understand are the rapid rise followed by slow approach to saturation of the signal and the width of the resonance lineshape (in the detuning Δ = ϵ′_{p} + ϵ_{p} − ϵ′_{s} − ϵ_{s}) as a function of time. The problem is not trivial because an initial ss′ → pp′ transition is followed by the propagation of sp → ps and s′p′ → p′s′ fluctuations, which are responsible for the broadening of the resonance. Although we address specifically the Rb experiments, the theory we present here is broadly applicable to analogous processes in any frozen gas of excited atoms.

In an earlier paper [5] we developed the general formalism, obtained exact results for the case when resonance broadening is negligible, and introduced simple approximations for the broadening. The same formalism and notation are used in this paper. However, in Ref. [5] we introduced the down and up states of an effective spin to describe the states s and p, and similarly for s′ and p′, but in this paper we prefer to describe an s → p transition as the creation of a Frenkel exciton [6–8], which then propagates through the medium. In this language, the initial creation of a pp′ pair launches two excitons, of types sp and s′p′, but we consider only the limit where the propagation of the s′p′ exciton is relatively unimportant and can be neglected. We prefer the exciton language because the literature on excitons in solids and liquids generally considers the dipole-dipole forces to be of major importance, as they are in a frozen Rydberg gas, and this literature has been useful to us in developing approximate analytical treatments. On the other hand, in spin systems the dominant role is usually played by the short-ranged exchange interaction and dipole-dipole forces are given comparatively little attention.

In this paper we present the results of numerical simulations for the case where N′ = 1 and N is large, which corresponds to the situation where a single s′ atom is initially present and a single sp exciton is launched by the resonant ss′ → pp′ transition. The calculated N_{p′} can be directly compared to the experimentally determined N_{p′}/N′ in the limit N′ ≪ N. When exciton propagation is suppressed, exact results are obtained by the methods of Ref. [5] and serve as a check on the numerical simulations.

We also develop analytical approximations for the average amplitude of a resonant transition and compare them with the results of the simulations. We find that the Cayley tree approximation agrees well with the simulations for a broad range of conditions, and is much better than the simple Lorentzian approximation introduced in Ref. [5]. The details of the numerical simulations are briefly outlined in Section II. Analytical approximations for the transition amplitude and for the exciton spectral density are developed in Section III. Remarkably, explicit formulas are obtained in the Cayley tree approximation. Representative graphs of the results, both numerical and analytical, are shown and discussed in Section IV, which is followed by the Conclusions.

II. DESCRIPTION OF NUMERICAL SIMULATIONS

Before discussing how the simulations are carried out, we review the basic elements of the theory presented in Ref. [5]. We consider one atom at the origin, initially in the state s′, in interaction with a gas of s atoms through a resonant ss′ → pp′ dipole-dipole interaction V. We also allow for interaction among the unprimed atoms via an sp → ps process mediated by another dipole-dipole potential U. This system is described by Eqs. (24a) and (24b) of Ref. [5], which are

i ȧ_{0} = Δa_{0} + ∑_{k=1}^{N} V_{k}c_{k},  (1a)

i ċ_{k} = V_{k}a_{0} + ∑_{l=1}^{N} U_{kl}c_{l}.  (1b)

Here a_{0}(t) is the amplitude of the state in which the atom at the origin is in state s′ and all other atoms are in state s, while c_{k}(t) is the amplitude of the state in which the atom at the origin is in state p′ and the atom at r_{k} is in state p, while all the others remain in state s. The quantities V_{k} and U_{kl} are the interaction potentials and Δ = ϵ′_{p} + ϵ_{p} − ϵ′_{s} − ϵ_{s} is the detuning from resonance. The potential V_{k} is of the dipole-dipole form

V_{k} = −\frac{2μμ′}{r_{k}^{3}} P_{2} (cos θ_{k}) ,  (2)

where θ_{k} is the angle between r_{k} and the magnetic field in the MOT used to confine the gas.^{1} Similarly U_{kl} is of the form

1The P_{2} (cos θ) angular dependence is also used, for instance, in Ref. [9]. In the absence of a magnetic field, each atom is in a multiplet of states and the angular dependences of each V_{k} and U_{kl} are described by matrices.

U_{kl} = −\frac{2μ^{2}}{r_{kl}^{3}} P_{2} (cos θ_{kl}) ,  (3)

for l ≠ k. We define U_{ll} = 0, so that the restriction l ≠ k is automatic in sums of the type appearing in Eq. (1b). Here μ is the dipole matrix element connecting the s and p states, and μ′ is the dipole matrix element connecting the s′ and p′ states. The atoms are assumed to be sufficiently cold that during the time scale of interest they move only a very small fraction of their separation, and therefore V_{k} and U_{kl} can be taken to be independent of time.

If we consider a column vector C (t) = (a_{0} (t) , c_{1} (t) , c_{2} (t) , …, c_{N} (t))^{T}, then clearly we can construct a matrix H such that the set of N +1 coupled differential equations (1) can be written in the form i Ċ = HC. If we let {λ_{m}} and {θ_{m}} be the eigenvalues and eigenvectors (represented as column vectors) of the real Hermitian matrix H, respectively, and if we further define e_{n} to be the column vector of length N +1 with the nth entry set equal to one and all others set equal to zero, then we can write e_{n} = ∑_{m=0}^{N} β_{nm}θ_{m} and θ_{m} = ∑_{n=0}^{N} β_{mn}e_{n}, where β_{nm} = θ_{m}^{T} · e_{n}. The system starts off in the state described by e_{0}, so it follows that

C (t) = exp (−iHt) e_{0}
= ∑_{m=0}^{N} β_{0m} exp (−iλ_{m}t) θ_{m}

= ∑_{m=0}^{N} ∑_{n=0}^{N} β_{0m}β_{mn} exp (−iλ_{m}t) e_{n}.  (4)

The quantity a_{0} (t) is just the e_{0} component of this expression, so

a_{0} (t) = e_{0}^{T} · C (t) = ∑_{m=0}^{N} β_{0m}β_{m0} exp (−iλ_{m}t) .  (5)

Thus we see that, in this form, solving for a_{0} at a given time becomes just a matter of diagonalizing the matrix H. To carry out the numerical simulations, we first generate a set of random positions for the unprimed atoms. For the results given in this paper these random positions are uniformly distributed, but it is just as easy to compute for a gas blob with a Gaussian density profile, as is the case in the experiments [1,2,4].^{2} From these positions the potentials V_{k} and U_{kl} and the matrix H are generated. We then diagonalize H and construct the quantity a_{0} at a set of values of t using Eq. (5). We do this for the same set of times for many random distributions of the unprimed atoms and thereby obtain an average of the quantities a_{0} and |a_{0}|^{2} at a given set of times over many configurations of the system. The signal measured by the experiments is then proportional to the average of 1 − |a_{0}|^{2}. One can directly generate Fourier transforms from Eq. (5), as we do in Section IV. We can also obtain plots of the linewidth as a function of time by generating signal versus time curves for many values of the detuning, Δ, and then determining for each time value which of these signals is equal to one half the value of the resonance signal.

2Results that directly compare calculations and experiments will be reported separately.

III. DESCRIPTION OF ANALYTICAL APPROXIMATION

A. Justification

In developing approximations, it is useful to cast the problem in the standard Green’s function language [10], as discussed in Section 3 of Ref. [11]. The Green’s function of Eq. (1b) satisfies the equation

ωG_{ln} = δ_{ln} + ∑_{m=1}^{N} U_{lm}G_{mn},  (6)

and then

a_{0} (ω) = \frac{i}{ω − Δ − ∑_{lm} V_{l}G_{lm}V_{m}} .  (7)

Actually, −ia_{0} is the 00 element of the Green’s function for the entire set of N +1 equations. Analogously, then, the on-site Green’s function for Eq. (1b) can be written in the form

G_{nn} = \frac{1}{ω − ∑_{lm} U_{nl}G_{lm,[n]}U_{mn}} ,  (8)

where G_{lm,[n]} is the Green’s function for a medium in which the n atom is absent. It is important to keep in mind that these Green’s function equations are exact.

The Cayley approximation [12,13] keeps only the diagonal (l = m) terms in the sum over l and m in Eqs. (7) and (8). Further, it assumes that each G_{ll,[n]} is independently distributed, and that there are sufficiently many sites that G_{ll,[n]} and G_{ll} have the same distribution.

We choose to set ω = iα and work in Laplace space instead of Fourier space because this makes many of the integrals that follow more obviously convergent. Up to a point, the same analysis can be carried out with the Fourier analogs, but Laplace transforms are needed to solve the Cayley equations exactly in part B of this Section. We also work with the inverses of these Green’s functions, setting a_{0} (iα) = 1/f_{0} (α) and G_{kk} (iα) = i/f_{k} (α), so in the Cayley approximation we have

f_{0} = α + iΔ + ∑_{k=1}^{N} \frac{V_{0k}^{2}}{f_{k}} ,  (9a)

f_{k} = α + ∑_{l=1}^{N} \frac{U_{kl}^{2}}{f_{l}} .  (9b)

B. U Process Only

First we consider just the U process and solve Eq. (9b), which is to say we look at the sp exciton by itself. Note that f_{k} is always real and positive. The distribution of f_{k} is given by

P (f_{k}) = ∫ ⟨δ (f_{k} − α − ∑_{l=1,}^{N} \frac{U_{kl}^{2}}{f_{l}} )⟩ ∏_{m≠k,m=1}^{N} P (f_{m}) df_{m},  (10)

where for a uniform gas of volume Ω,

⟨X⟩ = \frac{1}{Ω^{N}} ∫ ∏_{j=1}^{N} d^{3}r_{j} X.  (11)

Writing the delta function as a Fourier transform, we see that

P (f_{k}) = ∫_{−∞}^{∞} \frac{dq}{2π} ∫ ∏_{m≠k,m=1}^{N} P (f_{m}) df_{m} e^{iq(f_{k}−α)} ⟨exp (−iq ∑_{l=1}^{N} \frac{U_{kl}^{2}}{f_{l}} )⟩

= ∫_{−∞}^{∞} \frac{dq}{2π} e^{iq(f_{k}−α)} ∫ ∏_{m≠k,m=1}^{N} P (f_{m}) df_{m} \frac{1}{Ω^{N−1}} ∫ ∏_{l≠k,l=1}^{N} d^{3} r_{l} exp (−iq \frac{U_{kl}^{2}}{f_{l}} ) . (12)

Using the fact that (∏_{m} P_{m}) (∏_{l} E_{l}) = ∏_{m} P_{m}E_{m}, we can write

P (f_{k}) = ∫_{−∞}^{∞} \frac{dq}{2π} e^{iq(f_{k}−α)} ∏_{m≠k,m=1}^{N} ∫ P (f_{m}) df_{m} \frac{1}{Ω} ∫ d^{3} r_{m} exp (−iq \frac{U_{km}^{2}}{f_{m}} )

= ∫_{−∞}^{∞} \frac{dq}{2π} e^{iq(f_{k}−α)} [∫_{0}^{∞} df_{1} P (f_{1}) \frac{1}{Ω} ∫ d^{3}r_{1} exp (−iq \frac{U_{k1}^{2}}{f_{1}} )]^{N−1}

= ∫_{−∞}^{∞} \frac{dq}{2π} e^{iq(f_{k}−α)} (1 − ∫_{0}^{∞} df P (f) \frac{1}{Ω} ∫ d^{3}r \{1 − exp [−iq \frac{U^{2} (r)}{f}]\})^{N−1} , (13)

where in the last step we have simply changed variables from r_{1} to r = r_{1} − r_{k}. Using the averaging methods detailed in Eqs. (8–10) and the Appendix of Ref. [5], we see that

P (f) = ∫_{−∞}^{∞} \frac{dq}{2π} e^{iq(f−α)} exp [−u\sqrt{iq}Q (α)] ,  (14)

where

Q (α) = ∫_{0}^{∞} \frac{df}{\sqrt{f}} P (f) ,  (15)

and

u = \frac{16π^{3/2}}{9\sqrt{3}} μ^{2} \frac{N}{Ω} .  (16)

Using

∫_{−∞}^{∞} \frac{dq}{2π} e^{iqA−\sqrt{iq}B} = \frac{1}{2\sqrt{π}} \frac{B}{A^{3/2}} exp (−\frac{B^{2}}{4A}) θ (A) ,  (17)

for A and B real, we have

P (f) = \frac{1}{2\sqrt{π}} \frac{uQ (α)}{(f − α)^{3/2}} exp [−\frac{u^{2}Q^{2} (α)}{4 (f − α)}] θ (f − α) .  (18)

and so

Q (α) = ∫_{α}^{∞} \frac{df}{\sqrt{f}} \frac{1}{2\sqrt{π}} \frac{uQ (α)}{(f − α)^{3/2}} exp [−\frac{u^{2}Q^{2} (α)}{4 (f − α)}]
= \frac{1}{\sqrt{α}} exp [\frac{u^{2}Q^{2} (α)}{4α}] erfc [\frac{uQ (α)}{2\sqrt{α}}] .  (19)

We also have

⟨\frac{1}{f} ⟩ = ∫_{0}^{∞} \frac{df}{f} P (f)
= \frac{1}{α} − \frac{\sqrt{π}}{2} \frac{uQ (α)}{α^{3/2}} exp [\frac{u^{2}Q^{2} (α)}{4α}] erfc [\frac{uQ (α)}{2\sqrt{α}}] ,  (20)

or

⟨\frac{1}{f} ⟩ = \frac{1}{α} [1 − \frac{\sqrt{π}}{2} uQ^{2} (α)] .  (21)

Thus we can solve numerically the transcendental relation for Q (α) given in Eq. (19) and plug it into Eq. (21) to find ⟨1/f⟩. It is important to note that Eq. (19) defines Q as an analytic function of α in the half-plane where Re α > 0. Thus ⟨1/f⟩ can be analytically continued, in the variable ω = iα, from the positive imaginary axis to the entire upper half-plane.

C. U and V Processes in the Absence of Δ

Now we allow the V process to be present, but we still ignore the detuning, Δ. We start with the equations

f_{0} = α + ∑_{k=1}^{N} \frac{V_{0k}^{2}}{f_{k}} ,  (22a)

f_{k} = α + ∑_{l=1,l≠k}^{N} \frac{U_{kl}^{2}}{f_{l}} .  (22b)

Repeating the analysis we did before, we have

P_{0} (f_{0}) = ∫ ⟨δ (f_{0} − α − ∑_{k=1}^{N} \frac{V_{0k}^{2}}{f_{k}} )⟩ ∏_{l=1}^{N} P (f_{l}) df_{l},  (23)

which leads to the form

P_{0} (f_{0}) = ∫_{−∞}^{∞} \frac{dq}{2π} e^{iq(f_{0}−α)} \{∫_{0}^{∞} df P (f) \frac{1}{Ω} ∫ d^{3}r exp [−iq \frac{V^{2} (r)}{f}]\}^{N} .  (24)

At this point it is clear from our previous work that this equation leads to

P_{0} (f_{0}) = \frac{1}{2\sqrt{π}} \frac{vQ (α)}{(f_{0} − α)^{3/2}} exp [−\frac{v^{2}Q^{2} (α)}{4 (f_{0} − α)}] θ (f_{0} − α) ,  (25)

where

v = \frac{16π^{3/2}}{9\sqrt{3}} μμ′ \frac{N}{Ω} ,  (26)

and Q (α) is still given by Eq. (19), and depends on u. Therefore we obtain the result

⟨\frac{1}{f_{0}} ⟩ = ∫_{0}^{∞} \frac{df_{0}}{f_{0}} P_{0} (f_{0})
= \frac{1}{α} − \frac{\sqrt{π}}{2} \frac{vQ (α)}{α^{3/2}} exp [\frac{v^{2}Q^{2} (α)}{4α}] erfc [\frac{vQ (α)}{2\sqrt{α}}] ,  (27)

which again can be analytically continued.

D. U and V Processes in the Presence of Δ

If we further consider the effect of a detuning, Δ, then we have the full set of Eqs. (9). The equation for f_{k} is unchanged from the previous case, but because the equation for f_{0} is now complex we have to consider a distribution in the two variables f_{0}^{R} = Re f_{0} and f_{0}^{I} = Im f_{0}. We have

P_{0} (f_{0}^{R}, f_{0}^{I}) = ∫ ⟨δ (f_{0}^{R} − α − ∑_{k=1}^{N} \frac{V_{0k}^{2}}{f_{k}} ) δ (f_{0}^{I} − Δ)⟩ ∏_{l=1}^{N} P (f_{l}) df_{l}.  (28)

Clearly the delta function involving f_{0}^{I} does not participate in the averaging or in the integrations over f_{l}. Therefore we can treat the part involving f_{0}^{R} just as we did before, and the f_{0}^{I} part will just be carried along. Hence we find

P_{0} (f_{0}^{R}, f_{0}^{I}) = \frac{1}{2\sqrt{π}} \frac{vQ (α)}{(f_{0}^{R} − α)^{3/2}} exp [−\frac{v^{2}Q^{2} (α)}{4 (f_{0}^{R} − α)}] θ (f_{0}^{R} − α) δ (f_{0}^{I} − Δ) .  (29)

It follows, then, that

⟨\frac{1}{f_{0}} ⟩ = ∫_{0}^{∞} df_{0}^{R} ∫_{−∞}^{∞} df_{0}^{I} \frac{1}{f_{0}^{R} + if_{0}^{I}} P_{0} (f_{0}^{R}, f_{0}^{I})
= \frac{1}{α + iΔ} − \frac{\sqrt{π}}{2} \frac{vQ (α)}{(α + iΔ)^{3/2}} exp [\frac{v^{2}Q^{2} (α)}{4 (α + iΔ)}] erfc [\frac{vQ (α)}{2\sqrt{α + iΔ}}] .  (30)

E. Other Approximations

It is also possible to make approximations other than the one we have chosen. For example, instead of working with Eq. (9b), one can instead choose to work with

f_{k} = α + ∑_{l=1,}^{N} \frac{U_{kl}^{2}}{f_{k}} ,  (31)

which gives explicitly

f_{k} = \frac{1}{2} (α + \sqrt{α^{2} + 4 ∑_{l} U_{kl}^{2}}) .  (32)

The standard averaging procedure then leads to

⟨\frac{1}{f} ⟩ = \frac{2}{u\sqrt{π}} + \frac{α}{u^{2}} [exp (\frac{u^{2}}{α^{2}} ) erfc (\frac{u}{α}) − 1] .  (33)

Alternatively, we can replace 1/f_{l} by its average, obtaining

f_{k} = α + ∑_{l=1}^{N} U_{kl}^{2} ⟨\frac{1}{f} ⟩ ,  (34)

which leads to

α ⟨\frac{1}{f} ⟩ + \frac{\sqrt{π}}{2} \frac{u}{\sqrt{α}}\sqrt{⟨\frac{1}{f} ⟩} exp (\frac{u^{2}}{4α} ⟨\frac{1}{f} ⟩) erfc (\frac{u}{2α}\sqrt{⟨\frac{1}{f} ⟩}) = 1.  (35)

Eqs. (31) and (34) correspond to a Hubbard band and a CPA approximation, respectively. We have examined these approximations and found that they give results that do not agree with our simulations as well as the results obtained by the Cayley approximation. This is not at all surprising, since in the Cayley equations (9) the on-site Green’s functions for the U problem are allowed to vary from site to site, whereas in Eqs. (31) and (34) they are not.

IV. RESULTS

In this section we present and discuss the results of the numerical simulations described in Section II in comparison, when possible, with exact analytical results or with approximations that are good enough to be used instead of the simulations. We require that the graph of an approximate result should not entail much more numerical labor than the corresponding graph from the numerical simulations. Further, we omit comparisons of Laplace transforms, because two similar-looking Laplace transforms can correspond to distinctly different results in the time domain and to different spectral densities in the (real) frequency domain. We found that only the Cayley tree approximation for ⟨a_{0}⟩, as developed in Section III, comes close to meeting all these criteria and can be considered as “good” for a broad range of conditions, although in special limits simpler approximations suffice. We do not have, at this time, a good approximation for ⟨|a_{0}|^{2}⟩.

All the analytical results, exact and approximate, are for a single primed atom in an ideal, infinite gas of unprimed atoms with purely dipolar interactions. The simulations presented here were carried out for N unprimed atoms randomly distributed in a cubic volume floating in vacuum, with a single primed atom at the center of the cube. One can change the position of the unprimed atom (to study surface effects), the shape of the volume, the boundary conditions, and the interaction. One can also consider a non-uniform gas blob, such as one with a Gaussian profile, and put in some correlation between the positions of the atoms. Here we focus on the uniform, ideal gas as a standard.

For compactness of presentation, all energy and inverse time variables are scaled by the effective strength of the ss′ → pp′ transition, v, defined in Eq. (26). We have then results as a function of time or frequency and of two parameters: Δ/v, where Δ is the detuning from resonance, and u/v, where u is the equivalent of v for the sp → ps transition and is defined in Eq. (16). We display results for the following values of u: u=0

If the potential U is set equal to zero in Eq. (1b), then Eqs. (1) can be solved and averaged analytically as is done in Ref. [5]. This case therefore provides a convenient way to check to see if our numerical routines are working properly. It also provides a nice way to get a feeling for how many unprimed atoms need to be included and over how many realizations of the system we need to average to obtain good statistical results. In this case, of course, the sp exciton does not travel. u=v/4

We report some results for this case in order to show the gradual transition from the non-propagating to the propagating exciton, and also because this value of u/v corresponds to the Rb case when there are, initially, a few s atoms in a gas of s′ atoms (N ≪ N′). In effect, one interchanges the roles of primed and unprimed atoms with respect to the case u = 4v (see below), so that the s′p′ exciton is now the one that propagates. u=v

Resonant coupling to the excitons is nearly optimal in this case. In addition, for Δ = 0 finding a_{0} is equivalent to finding the on-site Green’s function G_{00} for a gas of N + 1 atoms. Then − (1/π) Im G_{00} (ω) = (1/π) Re a_{0} (ω) is the spectral density for exciton propagation. This and other quantities related to exciton propagation can then be studied independently from the mechanism of resonant excitation. u=4v

This value is chosen because it is close to the value u = 3.96v that pertains to the Rb experiments [1,2,4]. It also serves to illustrate, in general, what happens when exciton propagation is a large effect.

A. Results for N_{p′}/N′

The four parts of Fig. 1 display the numerical simulations of 1 − ⟨|a_{0}|^{2}⟩, corresponding to the experimental signal N_{p′}/N′, as a function of vt for u/v = 0, 1/4, 1 and 4, and enough values of the detuning Δ to give an idea of the resonance width in each case. The FWHM at large t can be directly deduced from the intercepts of the graphs with the right edge of the figure, and is discussed more fully in subsection D.

In addition, Fig. 1a shows the exact results for u = 0 that were reported in Ref. [5]. Comparison with the numerical simulations, which were carried out for N = 100, shows that the infinite system is already well simulated, but that noticeable wiggles remain after averaging over 10000 realizations. The reason for this persistence of fluctuations is that most realizations do not look at all like the average. In fact, the variance with realization (not shown) is comparable to the average and in particular for small t the average is given by (\sqrt{π}/2) vt, while each realization is quadratic in t, as already discussed in Ref. [5]. It can be seen from the graphs that the linear rise law, (\sqrt{π}/2) vt, is independent both of Δ and of u and holds for a substantial range of vt.

B. Results for ⟨a_{0} (t)⟩

The four parts of Figs. 2 and 3 display the numerical simulations of Re ⟨a_{0}⟩ and Im ⟨a_{0}⟩ , respectively, as a function of vt for u/v = 0, 1/4, 1, and 4. Recall that ⟨a_{0}⟩ is the probability amplitude that resonant exciton creation has not occurred during the time t. The striking feature of these graphs is that, off-resonance, oscillations with an approximate period Δ persist for a long time, not only when the exciton does not propagate (u = 0), but even for u = 4v.

The graph for u = 0 can be obtained analytically by the same methods used in Ref. [5] for computing 1 − ⟨|a_{0}|^{2}⟩. We give here only the following result for the Laplace transform:

⟨a_{0} (iα)⟩ = \frac{1}{α + iΔ} − \frac{\sqrt{π}}{2} \frac{v}{(α + iΔ)^{3/2} \sqrt{α}} exp [\frac{v^{2}}{4α (α + iΔ)}] erfc [\frac{v}{2\sqrt{α (α + iΔ)}}] . (36)

From this, ⟨a_{0} (t)⟩ has been obtained by expanding in powers of 1/α and inverting the Laplace transform term by term. As in Ref. [5], various formulas involving special functions can also be obtained. For u = 0, the exact analytical results (not shown here) agree with the numerical simulations about as well as in Fig. 1a. The initial dependence on time is given by

⟨a_{0}⟩ ≃ 1 − \frac{\sqrt{π}}{2} vt − iΔt  (37)

Looking at the graphs, it seems that the initial decrease of Re ⟨a_{0}⟩ depends on u, but this perception only indicates that the range of validity of the linear approximation (37) is small. On the other hand, the graphs of ⟨|a_{0}|^{2}⟩ in Fig. 1 look linear and independent of u over a wide range of vt. Note also that |⟨a_{0}⟩|^{2} ≃ 1 − \sqrt{π}vt decreases initially twice as fast as ⟨|a_{0}|^{2}⟩ ≃ 1 − (\sqrt{π}/2)vt. Of course, it must be true that the (mod) square of the average is smaller than or equal to the average of the square, but the above result shows that |⟨a_{0}⟩|^{2} is not even a fair approximation to ⟨|a_{0}|^{2}⟩.

C. Results for the Spectral Density

Figs. 4 to 7 present graphs of Re ⟨a_{0} (ω)⟩ that correspond to the graphs of ⟨a_{0} (t)⟩ in Figs. 2 and 3. For u = v and Δ = 0, Re ⟨a_{0} (ω)⟩ is the exciton spectral density (times π); more generally, it is an excitation spectral density for the resonant process ss′ → pp′. Each Figure gives a comparison of the numerical simulation with the Cayley approximation, which in particular for u = 0 reduces correctly to the exact result of Eq. (36), obtained from Eqs. (27) and (19) by analytic continuation. We compare spectral densities, rather than time graphs, in part because obtaining ⟨a_{0} (t)⟩ for the Cayley approximation is laborious, and in part because the comparison of spectral densities is instructive, as we now discuss.

In the notation of Section II, the spectral density is obtained from Eq. (5) as

Re a_{0} (ω + iε) = Im ∑_{m=0}^{N} \frac{β_{0m}β_{m0}}{ω − λ_{m} + iε}  (38)

in the limits ε → 0 and N → ∞. In practice it is difficult to go beyond N = 1000, and one must choose ε to be fairly large to obtain a smooth plot for the spectral density where the eigenvalues are sparse, even after averaging over 10000 realizations. On the other hand, the distribution of eigenvalues for a random dipolar gas is highly structured near ω = 0, and if ε is too large this structure will be missed. Rather than spending much time to go to large N and refine the numerical technique, we set N = 100, plot the results for ε/v = 0.1 and 0.01, and let the eye extrapolate to the correct limit. It is remarkable that the Cayley approximation agrees best with the smallest value of ε/v, except for wiggles that are usually in regions of sparse eigenvalues. This is exactly how the exact spectral density is expected to compare with a calculation for finite N. A finite ε can also be viewed as a simulation of the decoherence introduced by atomic motions, or of losses to other channels, such as the decay of the excited atomic states due to blackbody radiation [14]. Thus, ε is almost equivalent to the γ of the Lorentzian broadening model of Ref. [5], the difference being that in Ref. [5] iγ is added to ω for the Fourier transform of Eq. (1b) only, while here iε is added to every ω. As ε increases, the highly structured spectral densities of the ideal frozen gas become increasingly similar to Lorentzian lines. Because of unavoidable losses and incoherences, the experimental spectral densities will resemble the results of the simulations for some finite ε.

For u = 0 (Fig. 4) the spectral density is quite different from a Lorentzian. At Δ = 0 there is a quasi-gap around ω = 0. As Δ increases the gap extends from ω = 0 to ω = Δ.

For u = v/4 (Fig. 5) the sharp features of the u = 0 spectrum are rounded off and the gap is beginning to be filled, but is still recognizable.

For u = v (Fig. 6) the spectrum at Δ = 0 has the appearance of a single, asymmetric line. At Δ = v and Δ = 2v a shoulder on the left of the line is the remnant of the gap of Fig. 4.

For u = 4v (Fig. 7) there is a single line for all values of Δ, which becomes increasingly sharp as u increases. The lineshape is not Lorentzian, but even a small additional broadening will make it nearly so, as indicated by the calculation for ε/v = 0.01. This sharp line leads to the damped sinusoidal oscillations seen in the u = 4v graphs of Figs. 2 and 3. The general phenomenon of line narrowing can be discussed using the simple Lorentzian model for the exciton band, in which Eq. (1b) of Ref. [5] is replaced by

iċ_{k} = V_{k}a_{0} − iγc_{k}.  (39)

Solving the equations in this limit, one finds a pole at

ω = −i\frac{γ ∑_{k} V_{k}^{2}}{γ^{2} + Δ^{2}} + Δ,  (40)

which gives a decreasing width as u increases, since γ is proportional to u. One cannot simply replace ∑_{k} V_{k}^{2} by its average, which does not even exist. Instead, one can Fourier transform the average to find the time dependence

exp (−iΔt − \sqrt{γ_{eq}t}) ,  (41)

where

γ_{eq} = \frac{v^{2}γ}{γ^{2} + Δ^{2}} ,  (42)

and v has been defined in Eq. (26). A similar argument was applied in Ref. [5] to find the time dependence of ⟨|a_{0}|^{2}⟩.

From the comparison with the simulations, it is apparent that a shortcoming of the Cayley approximation is that it gives a symmetric spectral density at resonance (Δ = 0), while the correct spectral density is slightly asymmetric, except for u = 0. Correspondingly, ⟨a_{0} (t)⟩ is not purely real, as can be seen in Fig. 3. For Δ ≠ 0, the spectral density at ω < Δ is underestimated by the Cayley approximation.

Clearly, the asymmetry of the spectral density at Δ = 0 comes from the terms with l ≠ m in the denominator of Eq. (7), which are neglected in the Cayley approximation. It is surprising that these terms, which correspond to exciton propagation other than back and forth on the Cayley tree, have such small effect on ⟨a_{0}⟩ at Δ = 0. At finite Δ, the underestimate of Re ⟨a_{0} (ω)⟩ for ω < Δ is probably due to the neglect of these same terms. It seems complicated to include these terms, even approximately, in the Cayley tree formulation, but we have already shown how this can be done in the “CPA” approximation [11].

D. Results for the Resonance Linewidth

The time-dependent resonance lineshape is obtained from the ordinates of the curves in Fig. 1 at a given vt. In general the lineshape is not Lorentzian, but we never find a split line with a minimum at resonance as reported in Ref. [3]. Thus it is fair to characterize the line by its FWHM. This cannot be directly compared with an experimental width, which is usually the FWHM of a Lorentzian fitted to a noisy experimental line, but it is indicative of trends. For a more accurate comparison with experiment, one should of course convolute the calculated lineshape with the broadening from other effects.

Fig. 8 shows one half of the FWHM, w, as a function of vt for u = 4v. One sees the typical decrease of w from a universal t^{−1} behavior at small vt towards an asymptotic value, which is almost reached at vt = 15. The t^{−1} law, which results simply from the uncertainty principle, has already been discussed in Ref. [5], where on the basis of a small Δ expansion it was estimated that w = \sqrt{12}/t for small t.

The near-asymptotic value of w at vt = 15 is shown in Fig. 9 as a function of u/v. Overall, the graph is roughly linear except at small u/v, although a careful examination shows that there are really two linear regimes with somewhat different slopes that meet at u/v ≈ 3.5. There is an unexpected minimum around u = v/4, which is not a numerical artifact, but can be seen directly from Fig. 1 by comparing the plot for u = v/4 with those for u = 0 and u = v. What happens is that the resonance line always grows taller and wider as u increases, but for small u its height (the signal at Δ = 0) grows faster than its width.

V. CONCLUSIONS

In this paper, we report results of ongoing work on the general question: what are the electrodynamics and the optics of a gas of excited atoms interacting via dipolar fluctuations, for times short compared to the atom-atom collision time? We are considering in particular the time evolution of a resonant ss′ → pp′ process when N atoms in state s and N′ in state s′ are initially present, and all atoms are randomly distributed. At least one of these two initial states must be an excited state for the process to occur resonantly. Other resonant processes can be mapped on this prototype. For instance, the process pp → ss′ discussed in Refs. [15] and [3] is handled by setting s′ = s and reversing the roles of s and p. If the atoms are in Rydberg states then they have known energies and dipole moment matrix elements and the dipole-dipole interaction is dominant, so there are no adjustable parameters in the theory. More importantly for the experiments, with a frozen gas of Rydberg atoms one can easily do time-resolved spectroscopy because the evolution time scale is measured in microseconds. However, for the theory it does not matter that the atoms are of Rb or Cs, since other atoms in excited states would behave similarly on the corresponding time scale. In the language of the solid-state and chemical physics literatures, our problem can be described as resonant energy transfer to Frenkel excitons.

We have carried out detailed and extensive calculations, truly “computer experiments”, for a resonant ss′ → pp′ process when a single s′ atom is initially present in a frozen gas of s atoms. Neglecting the propagation of the s′p′ exciton to keep things simple, we consider in effect the resonant excitation and propagation of an sp exciton. The results apply, in practice, whenever only a few excitons of a given type are created, so that exciton-exciton interactions are negligible.

We can compute anything we wish about this system, but in this paper we report representative results for two quantities only: the average probability amplitude that no exciton has been created, ⟨a_{0}⟩, and the average probability that an exciton has been created, 1 − ⟨|a_{0}|^{2}⟩. One striking result is that the initial transition rate for this process is independent of exciton propagation and of detuning from resonance for surprisingly long times, and can thus be explicitly evaluated. Another notable feature is that the transition amplitude at detuning Δ shows persistent oscillations with period Δ, in particular when exciton propagation is a strong effect, as for u = 4v. Correspondingly, the spectral density for exciton creation is strongly peaked at Δ, and becomes narrower as u increases. Thus, the universal phenomenon of line narrowing actually simplifies the problem in the experimentally interesting limit of strong coupling (u ≫ v).

With a wealth of numerical results at our disposal, we have attempted to develop approximate formulas that are not simply fits to a particular numerical result. We have obtained good approximations for ⟨a_{0}⟩, but not yet for ⟨|a_{0}|^{2}⟩, which is notoriously more difficult to obtain than ⟨a_{0}⟩ and quite different from |⟨a_{0}⟩|^{2}. To develop approximation schemes we have drawn, and will continue to draw, from the literature on wave functions in random systems, and in particular on earlier work on excitons in disordered solids and liquids [13]. The best results, in our case, are given by the Cayley tree approximation.

Most of the work on random systems deals with on-site disorder or short-range interactions. On the contrary, the dipolar interaction U (r) ∼ 1/r^{3} is long range and has the distribution P (U) ∼ 1/U^{2} for a random distribution of r. The distribution of the cumulative variable U, where U^{2} = ∑_{k} U^{2} (r_{k}), plays a major role in the theory. It is given by

P (U) = \frac{1}{\sqrt{π}} \frac{u}{U^{2}} exp (−\frac{u^{2}}{4U^{2}} ) ,  (43)

where u is given by the appropriate equivalent of Eq. (16).^{3} Note that P (U) is not Gaussian and the central limit theorem does not apply, because P (U) does not have an r.m.s. value and is not even normalizable. Therefore it is not surprising that, as shown by the simulations, the distribution of eigenvalues and other statistical properties of the matrix U_{kl} are very different from those of a standard random matrix. This implies that random system results obtainable by random matrix theory are not applicable to dipolar interactions in an ideal gas. Luckily, the distribution in Eq. (43) makes the Cayley tree equations for ⟨a_{0}⟩ exactly soluble in this case.

ACKNOWLEDGMENTS

Acknowledgment is made to the Thomas F. and Kate Miller Jeffress Memorial Trust for the support of this research. We also wish to thank Ziad Maassarani for many useful discussions.

3The angular dependence of U modifies the value of u, leaving unchanged the functional form of Eq. (43). See Ref. [5], Appendix A.

REFERENCES

[1] W. R. Anderson, J. R. Veale, and T. F. Gallagher, Phys. Rev. Lett. 80, 249 (1998).

[2] W. R. Anderson, Ph.D. Dissertation, University of Virginia, 1996 (unpublished).

[3] I. Mourachko, D. Comparat, F. de Tomasi, A. Fioretti, P. Nosbaum, V. M. Akulin and P. Pillet, Phys. Rev. Lett. 80, 253 (1998).

[4] J. R. Lowell, Ph.D. Dissertation, University of Virginia, 1998 (unpublished).

[5] J. S. Frasier, V. Celli, and T. Blum, Phys. Rev. A 59, 4358 (1999).

[6] F. Bassani, G. Pastori Parravicini, and R. A. Ballinger, Electronic States and Optical Transitions in Solids (Pergamon, New York, 1975). See in particular chapter 6.

[7] Robert S. Knox, Theory of Excitons (Academic, New York, 1963).

[8] D. L. Dexter and R. S. Knox, Excitons (John Wiley & Sons, New York, 1965).

[9] L. Santos, G. V. Shlyapnikov, P. Zoller, and M. Lewenstein, Phys. Rev. Lett. 85, 1791 (2000).

[10] E. N. Economou, Green’s Functions in Quantum Physics (Springer, Berlin, 1983).

[11] V. Celli and J. S. Frasier, preprint cond-mat/9908230.

[12] R. Abou-Chacra, P. W. Anderson, and D. J. Thouless, J. Phys. C 6, 1734 (1973).

[13] D. E. Logan and P. G. Wolynes, J. Chem. Phys. 87, 7199 (1987); D. E. Logan and P. G. Wolynes, Phys. Rev. B 36, 4135 (1987); D. E. Logan and P. G. Wolynes, Phys. Rev. B 31, 2437 (1985).

[14] T. F. Gallagher, Rydberg Atoms (Cambridge, New York, 1994). See in particular chapter 5.

[15] V. M. Akulin, F. de Tomasi, I. Mourachko and P. Pillet, Physica D 131, 125 (1999).

FIGURES 17 18 19 20 21

FIG. 1. Plots of 1 − ⟨|a_{0} (t)|^{2}⟩ for u = 0, u = 0.25v, u = v, and u = 4v for several values of the detuning Δ, as shown on each graph. For u = 0, the exact result is shown for comparison.

23 24 25 26

FIG. 2. Plots of Re ⟨a_{0} (t)⟩ for u = 0, u = 0.25v, u = v, and u = 4v for three values of the detuning Δ.

28 29 30 31

FIG. 3. Plots of Im ⟨a_{0} (t)⟩, for the same values of u and Δ as in Fig. 2.

32 33 34 35

FIG. 4. Plots of Re ⟨a_{0} (ω)⟩ for u = 0 and Δ = 0, Δ = v, and Δ = 2v. Exact results are compared with numerical simulations of Re ⟨a_{0} (ω + iε)⟩ for two values of ε.

37 38 39

FIG. 5. Plots of Re ⟨a_{0} (ω)⟩ for u = 0.25v and Δ = 0, Δ = v, and Δ = 2v.

40 41 42 43

FIG. 6. As in Fig. 5, but for u = v.

44 45 46 47

FIG. 7. As in Fig. 5, but for u = 4v.

48 49

FIG. 8. A plot of the width as a function of time for u = 4v.

50 51

FIG. 9. A plot of the width at vt = 15 as a function of the ratio u/v.

52