text:
# Glassy behaviour in a simple topological model

## 1 Introduction

Over the last few years, much work has taken place on the topic of the glass transition and the behaviour of supercooled liquids, both experimentally and through computer simulations. A great deal of this activity has involved testing the predictions of Mode-coupling Theory (MCT), which has proved thus far to reliably describe many key aspects of the dynamics of glass formers (for a review of these, see ). However, many questions remain, particularly with regard to the extent to which one might expect a certain system to satisfy the predictions of MCT.

The bulk of the numerical studies have taken place on binary Lennard-Jones models, hard sphere systems or amorphous silica, all of which involve a number of parameters that can be chosen to suit the study e.g. to emulate experimental systems, or to avoid crystallization etc. By contrast, the aim of this paper is to investigate a model with as few parameters as possible, in order to establish the extent to which the same features and behaviour can be reproduced. The model we use is particularly simple in that it is purely topological, involving no length scale at all. There is also the added advantage that it shows no tendency to crystallize. Whilst the model is the dual of a two-dimensional atomic model, the dynamical driving force is not the same as in that case, and there is no true phase transition.

The lack of spatial co-ordinates within the model invalidates the investigation of many familiar functions, such as the van Hove correlation function, the structure factor and the mean squared displacement. However, suitable alternatives have been identified which are, in keeping with the ethos of this paper, simpler to calculate.

## 2 The Model

The model we use is that introduced by Aste and Sherrington in : it is an amorphous two-dimensional tiling of cells, with 3 edges incident on a vertex and 2 cells incident on an edge. The average number of sides of the cells is constrained to be 6 by the Euler Theorem . These tilings are topologically stable.

If $`n_i`$ is the number of sides of cell $`i`$, the topological charge $`q_i=6n_i`$ measures the deviation from the hexagonal configuration i.e. $`q_i`$ is a measure of inhomogeneity. Topological charge is conserved under local rearrangements of the cells.

We define the energy as follows:

$$E=\underset{i=1}{\overset{N}{}}(6n_i)^2=\mu _2N$$
(1)
where $`N`$ is the total number of cells in the system. This quantity is naturally associated with the degree of inhomogeneity in the configuration: the ground state (a perfect hexagonal tiling) possesses zero energy.

The system is allowed to evolve through moves as illustrated in Figure 1, known in the foam literature as T1 moves. These are local topological rearrangements exchanging topological charge between four cells. Two adjacent cells lose a side each and move apart, and two second-nearest neighbours gain a side each and become nearest neighbours. The energy change associated with such a move on 4 cells $`v,w,x,y`$ (as shown in Figure 1) with sides $`n_v,n_w,n_x`$ and $`n_y`$ is:

$$\mathrm{\Delta }E(n_w,n_y;n_v,n_x)=2(2+n_w+n_yn_vn_x)$$
(2)
We use Glauber-Kawasaki dynamics, which allows evolution of the system even at zero temperature provided the move decreases the energy or leaves it unchanged. The probability $`P`$ of performing a T1 move is given by:

$$P(n_w,n_y;n_v,n_x)=\frac{1}{1+\mathrm{exp}\left(\beta \mathrm{\Delta }E(n_w,n_y;n_v,n_x)\right)}(1\delta _{n_v,3})(1\delta _{n_x,3})(1\delta _{w,y})$$
(3)
where $`\beta `$ is the inverse temperature. The first two $`\delta `$-functions forbid the production of two-sided cells, whilst the last $`\delta `$-function forbids the production of tadpoles (i.e. self-neighbouring cells); unless these forbidden formations are present at the start, they will never appear.

In the following simulations the system consists of $`N=9900`$ cells, with periodic boundary conditions. We exhibit results of simulations starting from configurations obtained by randomly performing $`10^4`$N T1 moves on a perfect hexagonal tiling. This is equivalent to running the system at $`\beta =0`$, and results in an extremely disordered network, with a value for $`\mu _2`$ of approximately 13. Time is measured in units of $`N`$ attempted moves.

## 3 Relaxation Dynamics

We study the temporal dependence of the energy of the system by quenching from a disordered network to $`T=2.0`$ and allowing the system to equilibrate at that temperature for $`10^4`$N attempted moves. Cooling is then carried out by waiting a time $`t=\gamma N`$ at each temperature decrement of $`\delta T=0.05`$.

The results (averaged over 3 runs) are shown in Figure 2. The solid line superimposed on the graph is the equilibrium curve: the computation of this is covered thoroughly in . We can clearly see characteristic glassy behaviour: the system fails to reach equilibrium even at slow cooling rates for temperatures less than $`T0.2`$, and displays the typical strong dependence of the energy on the cooling rate, as found at low temperatures in glasses.

We also study the correlations within the system. The choice of correlation function for this particular system is by no means obvious. In previous work on this topological model, a persistence function has been studied, which measures the fraction of the total cells which have NOT been involved in a T1 move . However, this provides only limited information about the system: by the very definition of the function, information about a cell is thrown away once it has been involved in a move. If the system was in a metastable state in which it performed then undid many T1 moves, this would not be revealed.

Instead we have have chosen to use self-correlation functions: the first of these is

$$C(\tau ,\tau +t)=\frac{_{i=1}^N\left(n_i(\tau )6\right)\left(n_i(\tau +t)6\right)}{_{i=1}^N\left(n_i(\tau )6\right)^2}$$
(4)
We add in passing that the results for this particular function differ negligibly in a qualitative sense from those for an energy-energy self-correlation function such as that used in conjunction with the backgammon model .

We shall first present the results in equilibrium, where the correlation functions are stationary and possess no $`\tau `$ dependence. Thus the correlator $`C(\tau ,\tau +t)`$ becomes $`C(t)`$, a function of $`t`$ only. In order to obtain stationarity, the system was allowed to equilibrate for $`60,000N`$ at each temperature before measurements were taken, and for $`200,000N`$ at $`\beta =4`$. The energy was found to have equilibrated in each case and no $`\tau `$ dependence was found. For $`\beta >4`$ it became impossible to achieve equilibrium within any practical timescale; even if the curves appeared stationary, the energy had not equilibrated.

Figure 3(a) shows the correlation functions for a range of temperatures. At low values of $`\beta `$ (high temperatures) the correlator decays directly to zero. As the value of $`\beta `$ is increased, we see a shoulder develops, indicating that two-step relaxation is taking place. If we were able to attain stationarity at higher values of beta we would see this shoulder broaden to a plateau. The plateau can be seen more clearly in Figure 3(b), which shows the correlation function for higher values of $`\beta `$, with the waiting time of $`60,000N`$; however, one must keep in mind that this figure shows non-equilibrium results. Figures 3(c) and 3(d) show the strong $`\tau `$ dependence that exists out of equilibrium.

We can study the correlation function within the stationary regime. The relaxation around the shoulder or plateau, and the early stages of the departure from it, are known as the $`\beta `$-relaxation regime. Mode-coupling Theory predicts the behaviour of conventional correlation functions in the late part of this regime (i.e. during the plateau and the departure from it) to be governed by a von Schweidler law. This states that in stationary conditions the correlator $`C(t)`$ can be described as follows:

$$C(t)=fB\left(t/t_r(T)\right)^b$$
(5)
where $`f`$ is the plateau height (also known as the non-ergodicity parameter), and both $`B`$ and $`b`$ are positive constants. All three should be temperature-independent, and in addition $`b`$ should be independent of the choice of correlator. However, this theory is based on atoms interacting via a two-body potential and studies conventional correlation functions with explicit position dependence, such as the density-density correlator (for a recent review see ). Thus it is interesting to plot our rather unusual correlation function against rescaled time $`t/t_r`$ to investigate the predictions of MCT in our case.

Following the procedure of Kob , we define the relaxation time $`t_r`$ in Equation 5 as the time at which $`C(t)`$ first drops below $`\mathrm{e}^1`$. Figure 4 shows the correlation functions for a range of beta against rescaled time. All the functions collapse onto a master curve in the region after the plateau or shoulder, thus obeying the time-temperature superposition principle of MCT. In the late $`\beta `$-relaxation regime this curve can be fitted by a von Schweidler law with $`f=0.86`$, $`B=0.483`$ and $`b=0.42`$ as shown on Figure 4. The data fits well to this law, although because of the noise the margin of error on the fitting parameters is large.

Figure 5(a) shows the relaxation time against inverse temperature. This can be fitted extremely well by the following offset Arrhenius function:

$$t_r=A+B\mathrm{e}^{C/T}$$
(6)
where $`A,B,C`$ are constants. This is plotted on Figure 5(a), with $`A=3.13,B=0.0195`$ and $`C=2.85`$. Thus at very low temperatures (high inverse temperature) the system is displaying Arrhenius behaviour, which is characteristic of strong glass formers.

In the high temperature regime, one can fit the Vogel-Fulcher law typical of fragile glass-formers i.e. $`t_r=A\mathrm{e}^{C/(TT_o)}`$ with $`T_o=0.202`$ (see Figure 5(b)), but it does not suit the data as well as the offset Arrhenius curve.

The relaxation time is also subject to predictions by MCT, which states it should behave as:

$$t_r(TT_c^{MCT})^\gamma $$
(7)
close to a critical Mode-coupling temperature $`T_c^{MCT}`$, with $`T_c^{MCT}`$ defined through this equation. By performing such a fit on our data, we find that $`T_c^{MCT}=0.218`$, as shown on Figure 5(c), with $`\gamma =2.86`$. Neither this law nor the Vogel-Fulcher law fit the data well very close to $`T_c^{MCT}`$, although they both fit reasonably well slightly above $`T_c^{MCT}`$. However, it would seem that whilst one can fit both of these functions in certain regions, the best description at all temperatures is provided by the offset Arrhenius law.

Mode-coupling Theory predicts a relationship between the exponent $`\gamma `$ in Equation 7 and the exponent $`b`$ in Equation 5. They are linked by:

$$\gamma =\frac{1}{2a}+\frac{1}{2b}$$
(8)
where:

$$\frac{\mathrm{\Gamma }(1+b)^2}{\mathrm{\Gamma }(1+2b)}=\frac{\mathrm{\Gamma }(1a)^2}{\mathrm{\Gamma }(12a)}$$
(9)
Using our value of $`b=0.42`$, we find a value for $`\gamma `$ of $`3.12`$, which can be compared with the power law fit of $`\gamma =2.86`$. These two results agree to within $`10\%`$, which is reasonable given that the noise of the data makes it is very difficult to fit either exponent accurately.

## 4 Response Functions

A clear indication of the presence of aging is the violation of the fluctuation-dissipation theorem , which we shall briefly review. In equilibrium the response $`R_A`$ of an observable A to an applied conjugate field relates to the appropriate time-correlation function $`C_A=<A(t)A(0)>`$ as follows (in units of $`k_B=1`$):

$$R_A(t)=T^1\frac{C_A(t)}{t}$$
(10)
Out of equilibrium the relationship must be generalised since the correlation function is now a two-time quantity:

$$R_A(\tau ,t+\tau )=T^1X_A(\tau ,t+\tau )\frac{C_A(\tau ,t+\tau )}{\tau }$$
(11)
with $`X_A(\tau ,t+\tau )`$ defined by the above equation. If $`X_A(\tau ,t+\tau )`$ is equal to 1, equation 11 reverts to the equilibrium case.

In order to investigate autocorrelations, one can apply a field at time $`\tau `$ which has a magnitude $`h_o`$, but is randomly positive or negative across the entire system. That field is left switched on, and the integrated response $`G(\tau ,t+\tau )`$ is measured, where:

$$G(\tau ,t+\tau )=h_o_\tau ^{t+\tau }R_A(t^{},t+\tau )𝑑t^{}$$
(12)
For many systems it has been observed that for $`t`$ and $`\tau `$ both large, $`X_A(\tau ,t+\tau )`$ depends on $`t,\tau `$ only through the correlator i.e. $`X_A(\tau ,t+\tau )=x\left(C(\tau ,t+\tau )\right)`$. If that is indeed the case, one finds:

$$\frac{TG(\tau ,t+\tau )}{h_o}=_1^Cx(C^{})𝑑C^{}$$
(13)
Under those conditions the slope of a parametric plot of $`TG(\tau ,t+\tau )/h_o`$ against the relevant correlation function is the interesting quantity, since:

$$\frac{}{C}\left(\frac{TG(\tau ,t+\tau )}{h_o}\right)=x(C)$$
(14)
One would expect to find a slope of exactly $`1`$ where the fluctuation-dissipation ratio is upheld. Where the ratio is broken, the slope gives information about the form of $`X_A(\tau ,t+\tau )`$.

In our particular case, we apply a perturbation $`h_o_{i=1}^Nϵ_i(n_i6)`$ at time $`\tau `$. $`ϵ_i`$ is randomly assigned to be $`+1`$ or $`1`$, and $`h_o`$ is chosen carefully to ensure linear response whilst also obtaining a reasonable signal-to-noise ratio. We follow the integrated response $`G(\tau ,t+\tau )=_{i=1}^Nϵ_i(n_i6)`$. It is easy to see that the appropriate correlation function for this observable is of the form of that given in equation 4, although one must be careful with normalisation factors. The data is in general very noisy and it is necessary to average over at least 8-10 different field configurations to be able to reach any conclusions.

Before presenting results, we draw attention to the unusual nature of our integrated response function. At low temperatures almost all the cells within the system are 5, 6 or 7-sided, with the vast majority being 6-sided. This means that a large proportion of the system makes no contribution to the response whatsoever. Furthermore, processes which reduce the energy by turning pentagons and heptagons into hexagons in fact reduce the response. This will be discussed in more detail in the following paragraphs; in the meantime, it is sufficient to remember that the integrated response function used here has some unusual features.

Figure 6 shows the results for $`\beta =4`$ and $`5`$. These are non-equilibrium results in all cases, as the waiting times are not long enough for the system to have equilibrated. The superimposed straight lines have slope -1, but are not best fits; they are there to give the eye a comparative reference. These plots should be examined from right to left i.e. short times are on the right, where $`C(\tau ,\tau +t)`$ is close to $`1`$, and long times are on the left. We see that for short times the fluctuation-dissipation ratio is upheld; the slope is -1 for all values of $`\tau `$, although the intercept does of course vary. In fact, the ratio is not broken until $`C(\tau ,\tau +t)`$ has decreased to $`0.6`$: this is considerably lower than the value of $`C(\tau ,\tau +t)`$ at which the plateau occurs. Thus the fluctuation-dissipation ratio is obeyed for a timescale larger than that of the onset of the plateau. After the ratio has been broken, curves for different $`\tau `$ behave differently, tending towards the equilibrium straight line of slope $`1`$ as $`\tau `$ increases. The non-monotonicity displayed clearly in Figure 6(a) occurs in every case, although at lower temperatures one cannot run long enough to see it, as in Figure 6(b). Non-monotonicity of the response has been noted in several other models and will be discussed briefly in the next section. In the meantime, with regards to this model this unusual behaviour can be understood if we now turn to a discussion of the processes dominating the evolution of the system.

At low temperatures, out of equilibrium, moves which increase the energy occur very rarely; we can consider the system to evolve by a combination of zero-energy moves and moves which reduce the energy. After a short time ($`10^3N`$) the system consists almost entirely of 5, 6 and 7-sided cells. Under these conditions the possible moves are as shown in Figure 7(a) and 7(b). The moves in Figure 7(c) are considered very unlikely due to the low probability of finding four non-hexagons clustered together; this is corroborated by the extreme rarity of finding octagons and rectangles present. We see from Figure 7(a) that in certain topological arrangements, pentagons and heptagons can effectively diffuse freely through the system; in particular, pentagons and heptagons can ‘pair up’ to move through the system together with no energy cost. To reduce the energy, one must annihilate two or more ‘defects’ (i.e. non-hexagons) as in Figure 7(b); this requires the absorption of a 5-7 pair by an isolated pentagon or hexagon, or the annihilation of two 5-7 pairs. This brings us to a conceptual picture containing both fast and slow dynamics: the fast dynamics is due to the rapid diffusion of 5-7 pairs moving freely through the system, whereas the slow dynamics is due to annihilation processes. The latter can be broken down into two different types: absorption of a 5-7 pair by either a heptagon or a pentagon (which can also rearrange the network such that isolated pentagons and heptagons become 5-7 pairs), and the complete annihilation of two 5-7 pairs. Of course, in equilibrium an equal number of 5-7 pairs are created as annihilated/absorbed, but in the relaxation to equilibrium, annihilation has the upper hand as the system starts from a configuration in which there are more 5-7 pairs present.

The effect of the field is to try and ‘pin’ defects onto the appropriate cell i.e pentagons try to settle on sites with $`ϵ_i=+1`$ and heptagons on $`ϵ_i=1`$. Thus with the field switched on, the diffusion processes move the defects around until they are on the appropriate sites, and in this way increase the value of the response $`TG(\tau ,\tau +t)/h_o`$. However, annihilation processes remove defects from the system entirely, thus reducing the response and also reducing $`C(\tau ,\tau +t)`$. We have competition between the annihilation processes and the diffusive processes; when annihilation dominates, the slope on the parametric plot becomes positive, as seen for small values of $`C(\tau ,\tau +t)`$ in Figure 6(a). If one looks instead at a quantity such as $`TG(\tau ,\tau +t)/h_o\mu _2`$ (in some sense a ‘response per defective cell’) one finds that the slope remains negative: this tells us that the diffusion processes are succeeding in placing more of the defects that remain on the appropriate sites, but the number of those defects present in the system is declining. The value of the response reached when the system finally achieves equilibrium depends on the waiting time $`\tau `$, and increases as $`\tau `$ increases; this is because both the response and the correlation function have been normalised by a factor that is equivalent to the energy at time $`\tau `$ (see equation 4).

This picture of diffusion and annihilation/creation processes is helpful to an understanding of the plateaux in the correlation functions at low temperatures. The initial descent to the plateau is due to diffusing 5-7 pairs, which quickly move the network away from the starting configuration. However, the lone defects are trapped, and can only be freed by a move which costs energy or by absorption of a 5-7 pair; both of these processes occur on a timescale which is temperature-dependent, and thus the length of the plateau itself is dependent on the temperature. This picture of trapped or caged defects is conceptually similar to that typically used when dealing with Lennard-Jones binary models .

## 5 Concluding remarks

In summary, we have studied a simple two-dimensional topological glassy model, testing some of the predictions of Mode-coupling theory for this system. We find that this theory can indeed be used to describe the behaviour of the system in so far as we have investigated, although this description fits some features more closely than others. The correlation function follow a von Schweidler law as predicted in the late $`\beta `$-relaxation regime, and the exponent of this roughly agrees with the exponent of the power law fitted to the relaxation time data. However, the relaxation times fit better to an offset Arrhenius law in all temperature regimes.

Investigation of the fluctuation-dissipation ratio reveals an interesting feature, namely non-monotonicity of the response. This is a feature which has also been observed in other simple models; for example, vibrated granular media , constrained Ising chains and the Backgammon model in one dimension , and also recently in a two-dimensional short-range spin model with uniform ferromagnetic 3-body interactions . In each case, the dynamics at low temperatures can be considered to involve activation over energy barriers; as mentioned in , this raises the question as to whether this is a generic feature of models with activated dynamics.

## 6 Acknowledgements

The authors would like to thank T. Aste and J. P. Garrahan for helpful discussions, and EPSRC(UK) for financial support - DS for research grant GR/M04426, and LD for research studentship 98311155.


format: text

text:
# Deeply Virtual Compton Scattering and Skewed Parton Distribution

## Introduction

The study of the structure of the nucleon is one of the most important frontiers in strong interaction physics. There are still many unanswered questions largely due to the non-perturbative nature of the bound state problem in QCD. Two traditional types of observables have been studied extensively, both theoretically and experimentally, for the last forty years: the parton (quark and gluon) distribution functions (PDF’s) (via deeply inelastic scattering (DIS) or Drell-Yan processes), and elastic form factors of the nucleon. In the past few years, studies of a new type of nucleon observable, the skewed parton distributions (SPD’s), have flourished(eg,spd ). The SPD’s generalize and interpolate between the ordinary PDF’s and elastic form factors and therefore contain rich structural information. They can be measured in (exclusive) diffractive processes in which the nucleon recoils elastically after receiving a non-zero momentum transfer in the so-called deeply virtual limit111$`Q^2\mathrm{}`$ while keeping Bjorken variable $`x=x_B=\frac{Q^2}{2pq}`$ fixed, and $`Q^2>>t=(p^{}p)^2`$ where $`p`$ ($`p^{}`$) is the initial(recoil) nucleon momentum..

## Theory: From DIS & PDF’s to DVCS & SPD’s

Ordinary PDF’s in DIS are accessed through an optical theorem that relates the imaginary part of a forward (Compton) scattering amplitude to the cross section. Through operator product expansion (OPE) one factorizes the hard scattering from the soft physics. The hard scattering can be calculated order by order in perturbation theory while the soft part is parametrized as PDF’s. They are essentially matrix elements of light-cone bilocal (quark and gluon) operators between equal momentum (symmetric) states. The factorization scale dependence of the matrix elements/PDF’s is governed by a renormalization group equation. The cross section is related to (usually a linear combination of) PDF’s.

The SPD’s are, on the other hand, essentially matrix elements of the same light-cone bilocal operators between different momentum (asymmetric) states. They can be accessed thus via deeply virtual processes like diffractive vector meson production ($`epep^{}VM`$)BrodskyVM ; Radyphlt962 ; Collins&FS97 ; Frankf&Strikman98  and deeply virtual Compton scattering (DVCS)Dmuel ; Ji97 ; Radyphlt961 . DVCS(See fig.1. (b) is the lowest order “handbag” diagram in an e-p collision, where $`x`$ and $`x^{}`$ denote the longitudinal momentum fractions of the interacting parton,) is a non-forward process signified by non-zero $`t`$ with the longitudinal momentum (fraction) transfer $`xx^{}\zeta =x_B`$ and also in general a non-zero transverse momentum transfer.

One still has valid factorization theoremCollins&FS97 ; Factorization  and OPEDmuel ; nfvcs  (fig.2) in the non-forward case. That is, the DVCS amplitude can be factorized into a hard scattering part calculable perturbatively (the upper part of the OPE diagram and also a crossed term) and a soft part (lower part of the diagram) parametrized as SPD’s.

SPD’s, like usual PDF’s, depend on factorization scale. Its QCD evolution is governed again by a renormalization group equation through a kernel that has been worked out to next-to-leading-order (NLO)Belitsky .222Also worked out is the NLO correction to hard scattering in DVCS. The QCD $`Q^2`$ evolution of the SPD’s has been studied extensively and has been shown to exhibit characteristics of both the DGLAP evolution of usual PDF’s and the ERBL evolution of meson distribution amplitudes, depending on different kinematic regions(eg, spd ; nfvcs ). In particular, in the small $`x_B`$, ie, small $`\xi =\frac{x_B}{2(1x_B/2)}`$ region SPD’s are closely connected with usual PDF’sCollins&FS97 . For $`x>>\xi `$, SPD’s $``$ PDF’s, while for $`x\xi `$, in leading $`\mathrm{log}\frac{1}{x}\mathrm{log}\frac{1}{\xi }`$, SPD’s $``$ forward PDF’sBrodskyVM  and at large factorization scale the $`Q^2`$ evolution tends to wash out effects of the asymmetry $`\xi `$ (eg, Martin ). The region $`x<\xi `$ is least well-known.

SPD’s are indeed form factors of the non-forward scattering amplitude $`\gamma ^{}p\gamma p^{}`$, eg, take quarks (similar for gluons)

$`FT`$ $`p^{},s^{}|\overline{\psi }(0)\gamma ^\mu \psi (z)|p,s=H(x,\xi ,t)\overline{u}(p^{},s^{})\gamma ^+u(p,s)`$
$`+E(x,\xi ,t)\overline{u}(p^{},s^{}){\displaystyle \frac{i\sigma ^{+\nu }(p^{}p)_\nu }{2M}}u(p,s)+\mathrm{},(\mathrm{Quark}\mathrm{Spin}\mathrm{Sum})`$
$`FT`$ $`p^{},s^{}|\overline{\psi }(0)\gamma ^+\gamma _5\psi (z)|p,s=\stackrel{~}{H}(x,\xi ,t)\overline{u}(p^{},s^{})\gamma ^+\gamma _5u(p,s)`$
$`+\stackrel{~}{E}(x,\xi ,t)\overline{u}(p^{},s^{}){\displaystyle \frac{(p^{}p)^+\gamma _5}{2M}}u(p,s)+\mathrm{},(\mathrm{Difference})`$
where $`H(x,\xi ,t)`$, $`E(x,\xi ,t)`$, $`\stackrel{~}{H}(x,\xi ,t)`$, $`\stackrel{~}{E}(x,\xi ,t)`$ are skewed quark distributions with $`2\xi =\frac{x_B}{1x_B/2}`$ 333Notations in the literature differ. That of Ji’s is used and mom. fractions refer to $`\frac{1}{2}(p+p^{})`$. and $`FT`$ denotes Fourier Transform.

While in general two SPDs correspond to one usual PDF, in the forward limit of $`p=p^{}`$ ($`\xi 0`$ and $`t0`$), SPD’s are reduced to the normal distributions: $`H(x,0,0)=q(x),\stackrel{~}{H}(x,0,0)=\mathrm{\Delta }q(x),`$ where $`q(x)`$ and $`\mathrm{\Delta }q(x)`$ are the conventional forward quark and quark helicity distributions (Similar equations hold for gluons). At the same time, the first moments of these SPD’s are related to nucleon form factors of corresponding EM or Axial currents by the sum rules:

$`{\displaystyle _1^1}𝑑xH(x,\xi ,t)=F_1(t),`$ $`{\displaystyle _1^1}𝑑xE(x,\xi ,t)=F_2(t),`$
$`{\displaystyle _1^1}𝑑x\stackrel{~}{H}(x,\xi ,t)=G_A(t),`$ $`{\displaystyle _1^1}𝑑x\stackrel{~}{E}(x,\xi ,t)=G_P(t).`$
At the same time, SPD’s have many new features, in contrast to the usual PDF’s.

* The non-forward amplitude to lowest $`O(\alpha _s)`$ is generally related to a integration over the SPD’s of type $`Amp𝑑x\frac{1}{x\zeta +iϵ}f(x,\zeta ,t)`$ ($`f`$ denotes SPD’s).
* Cross section is obtained by squaring the amplitude. There is not an optical theorem, nor a simple probability interpretation. SPD’s can be viewed as overlap of wavefunctions between different parton numbers (Fock states)Fock
* SPD’s are interference/correlation functions of different wave functions/ probability amplitudes. The extreme case with $`x>0`$ and $`x^{}<0`$ can be understood by re-interpret $`x^{}`$ line as antiparton with momentum fraction $`x^{}`$ resulting in “extracting $`q\overline{q}`$-pair” from the nucleon.

SPD(via DVCS) is the only place so far probing the orbital angular momentum of partonsJi97 . Moments of the SPD’s provide information on quark and gluon contributions to nucleon’s spin because they are closely related to the form factors of the (QCD) energy-momentum tensor. For example, one can measure the skewed quark distribution in spin-averaged experiments and extract form factors of the tensor. Extrapolating to $`t=0`$ one can obtain the total (spin+orbital) quark contribution to the nucleon spin $`J_q(0)`$ by taking $`t=0`$ in Ji’s sum rule:

$$\frac{1}{2}𝑑xx(H_q(x,\xi ,t)+E_q(x,\xi ,t))=J_q(t).$$

There is a similar sum rule for gluon and one finds $`J_q(0)+J_g(0)=\frac{1}{2}`$.

## Phenomenology and Experiments

Using known usual PDF’s as input, models of SPD’s have been studies by making ansatz that fulfills the (some of which very non-trivial) general properties of SPD’s(eg, spd ; models ). There are also studies on contributions from meson exchange in the ERBL region ($`x>0`$,$`x^{}<0`$, on SPD’s as overlap of proton wavefunctions for $`x>\xi `$Diehl  and on different dynamical models of SPD’s, including bag model, instanton vacuum and quark soliton model (eg, dmodels ).

Experiments usually access integrals of SPD’s with parton momenta that are multi-variable and can not be directly obtained from cross sections. Thus two main issues remain, one being experimental difficulty, the other extracting SPD’s from data.

The key background to DVCS is the QED Bethe-Heitler(BH) process(Compton scattering), which depend on EM form factor, and its interference with DVCS. At higher $`Q^2`$ the background is smaller, but so is the signal. This makes measuring DVCS cross section very difficult experimentally. The first experimental measurement come from HERAhera , where the preliminary on 96/97 data has shown evidence of the DVCS signal. However, the interference also makes possible exploring DVCS at the amplitude level, measuring its imaginary and real part independently.444A proposal for JLAB Hall Ajlab  proposed measuring cross section difference for leptons of opposite helicities, which is proportional to the interference of the imaginary part of the DVCS amplitude with a known BH weight. This quantity turns out to be a linear combination of SPD’s, thus if scaling is reached at the kinematic region with $`6Gev`$ beam, it would be a measurement of SPDs’ contribution. Also since the higher twist effects is only down $`\frac{1}{Q}`$ and has a different angular distribution (thus not masked by leading-twist) it can give an estimation of these effects. There has been a lot of data on vector meson production (eg, hermes ; heravm ), but the extraction of SPD’s is still a pending task.555Another reason of the interest in exclusive diffractive VM at small-x is that the cross section is predicted to be proportional to the square of gluon densityBrodskyVM . This could be a direct measurement of the glue in the proton, rather than getting it from evolution as in DIS. Table 1 is a very brief overview of the current experimental status.

Independent of the actual form of SPD’s, factorization predicts that in the scaling limit of $`Q^2\mathrm{};x_B,t`$ fixed, DVCS $`\frac{1}{Q^0}`$ and VM production $`\frac{1}{Q}`$. There are also helicity selection rules stating that for DVCS only $`\gamma ^{}(T)\gamma (T)\frac{1}{Q^0}`$ and $`\gamma ^{}(L)\gamma (T)`$ is power suppressed while for VM only $`\gamma ^{}(L)VM(L)\frac{1}{Q}`$ and all else are power suppressed. This is due to helicity conservation and introduces helicity parton distributions (eg, helicity flip gluon distribution that has no counterpart in DIS)helipdf , while being higher-twist gives the suppression. Therefore by looking at the angular distributions one have a window on higher-twist effects htwist  since it is not masked by leading twist because of different helicity.

Other SPD-related experimental processes been proposed include diffractive dijets (photo-production) and crossed DVCS (eg, $`\gamma ^{}\gamma \pi ^+\pi _{}`$)other , etc.

## Acknowledgements

I thank Xiangdong Ji and Macus Diehl for helpful discussions. This research is partly supported by a grant from Manhattanville College.


text:
# Resonant processes and exciton propagation in a frozen gas

## 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. and Mourachko et al.. The work presented here was motivated, in part, by the desire to understand those experiments and the subsequent work of Lowell et al. highlighting the dynamic aspects of these resonant, many particle systems. The measurements of Refs.  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 $`\mathrm{\Delta }=ϵ_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 $`spps`$ 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 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.  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 $`sp`$ transition as the creation of a Frenkel exciton, 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.  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. . 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. . 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 $`spps`$ process mediated by another dipole-dipole potential $`U`$. This system is described by Eqs. (24a) and (24b) of Ref. , which are

$`i\dot{a}_0`$ $`=`$ $`\mathrm{\Delta }a_0+{\displaystyle \underset{k=1}{\overset{N}{}}}V_kc_k,`$ (1a)
$`i\dot{c}_k`$ $`=`$ $`V_ka_0+{\displaystyle \underset{l=1}{\overset{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 $`𝐫_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 $`\mathrm{\Delta }=ϵ_p^{}+ϵ_pϵ_s^{}ϵ_s`$ is the detuning from resonance. The potential $`V_k`$ is of the dipole-dipole form

$$V_k=\frac{2\mu \mu ^{}}{r_k^3}P_2\left(\mathrm{cos}\theta _k\right),$$
(2)
where $`\theta _k`$ is the angle between $`𝐫_k`$ and the magnetic field in the MOT used to confine the gas.The $`P_2\left(\mathrm{cos}\theta \right)`$ angular dependence is also used, for instance, in Ref. . 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. Similarly $`U_{kl}`$ is of the form

$$U_{kl}=\frac{2\mu ^2}{r_{kl}^3}P_2\left(\mathrm{cos}\theta _{kl}\right),$$
(3)
for $`lk`$. We define $`U_{ll}=0`$, so that the restriction $`lk`$ is automatic in sums of the type appearing in Eq. (1b). Here $`\mu `$ is the dipole matrix element connecting the $`s`$ and $`p`$ states, and $`\mu ^{}`$ 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\left(t\right)=(a_0\left(t\right),c_1\left(t\right),c_2\left(t\right),\mathrm{},c_N\left(t\right))^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\dot{C}=HC`$. If we let $`\left\{\lambda _m\right\}`$ and $`\left\{\theta _m\right\}`$ 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 $`n`$th entry set equal to one and all others set equal to zero, then we can write $`e_n=_{m=0}^N\beta _{nm}\theta _m`$ and $`\theta _m=_{n=0}^N\beta _{mn}e_n`$, where $`\beta _{nm}=\theta _m^Te_n`$. The system starts off in the state described by $`e_0`$, so it follows that

$`C\left(t\right)`$ $`=`$ $`\mathrm{exp}\left(iHt\right)e_0`$ (4)
$`=`$ $`{\displaystyle \underset{m=0}{\overset{N}{}}}\beta _{0m}\mathrm{exp}\left(i\lambda _mt\right)\theta _m`$ (5)
$`=`$ $`{\displaystyle \underset{m=0}{\overset{N}{}}}{\displaystyle \underset{n=0}{\overset{N}{}}}\beta _{0m}\beta _{mn}\mathrm{exp}\left(i\lambda _mt\right)e_n.`$ (6)
The quantity $`a_0\left(t\right)`$ is just the $`e_0`$ component of this expression, so

$$a_0\left(t\right)=e_0^TC\left(t\right)=\underset{m=0}{\overset{N}{}}\beta _{0m}\beta _{m0}\mathrm{exp}\left(i\lambda _mt\right).$$
(7)
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.§§§Results that directly compare calculations and experiments will be reported separately. 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. (7). 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 $`\left|a_0\right|^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\left|a_0\right|^2`$. One can directly generate Fourier transforms from Eq. (7), 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, $`\mathrm{\Delta }`$, and then determining for each time value which of these signals is equal to one half the value of the resonance signal.

## III DESCRIPTION OF ANALYTICAL APPROXIMATION

### A Justification

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

$$\omega G_{ln}=\delta _{ln}+\underset{m=1}{\overset{N}{}}U_{lm}G_{mn},$$
(8)
and then

$$a_0\left(\omega \right)=\frac{i}{\omega \mathrm{\Delta }_{lm}V_lG_{lm}V_m}.$$
(9)
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}{\omega _{lm}U_{nl}G_{lm,[n]}U_{mn}},$$
(10)
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  keeps only the diagonal $`\left(l=m\right)`$ terms in the sum over $`l`$ and $`m`$ in Eqs. (9) and (10). 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 $`\omega =i\alpha `$ 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\left(i\alpha \right)=1/f_0\left(\alpha \right)`$ and $`G_{kk}\left(i\alpha \right)=i/f_k\left(\alpha \right)`$, so in the Cayley approximation we have

$`f_0`$ $`=`$ $`\alpha +i\mathrm{\Delta }+{\displaystyle \underset{k=1}{\overset{N}{}}}{\displaystyle \frac{V_{0k}^2}{f_k}},`$ (11a)
$`f_k`$ $`=`$ $`\alpha +{\displaystyle \underset{l=1}{\overset{N}{}}}{\displaystyle \frac{U_{kl}^2}{f_l}}.`$ (11b)
### 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\left(f_k\right)=\delta \left(f_k\alpha \underset{l=1,}{\overset{N}{}}\frac{U_{kl}^2}{f_l}\right)\underset{mk,m=1}{\overset{N}{}}P\left(f_m\right)df_m,$$
(12)
where for a uniform gas of volume $`\mathrm{\Omega }`$,

$$X=\frac{1}{\mathrm{\Omega }^N}\underset{j=1}{\overset{N}{}}d^3r_jX.$$
(13)
Writing the delta function as a Fourier transform, we see that

$`P\left(f_k\right)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dq}{2\pi }}{\displaystyle \underset{mk,m=1}{\overset{N}{}}P\left(f_m\right)df_me^{iq\left(f_k\alpha \right)}\mathrm{exp}\left(iq\underset{l=1}{\overset{N}{}}\frac{U_{kl}^2}{f_l}\right)}`$ (14)
$`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dq}{2\pi }}e^{iq\left(f_k\alpha \right)}{\displaystyle \underset{mk,m=1}{\overset{N}{}}P\left(f_m\right)df_m\frac{1}{\mathrm{\Omega }^{N1}}\underset{lk,l=1}{\overset{N}{}}d^3r_l\mathrm{exp}\left(iq\frac{U_{kl}^2}{f_l}\right)}.`$ (15)
Using the fact that $`\left(_mP_m\right)\left(_lE_l\right)=_mP_mE_m`$, we can write

$`P\left(f_k\right)`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dq}{2\pi }}e^{iq\left(f_k\alpha \right)}{\displaystyle \underset{mk,m=1}{\overset{N}{}}}{\displaystyle P\left(f_m\right)𝑑f_m\frac{1}{\mathrm{\Omega }}d^3r_m\mathrm{exp}\left(iq\frac{U_{km}^2}{f_m}\right)}`$ (16)
$`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dq}{2\pi }}e^{iq\left(f_k\alpha \right)}\left[{\displaystyle _0^{\mathrm{}}}𝑑f_1P\left(f_1\right){\displaystyle \frac{1}{\mathrm{\Omega }}}{\displaystyle d^3r_1\mathrm{exp}\left(iq\frac{U_{k1}^2}{f_1}\right)}\right]^{N1}`$ (17)
$`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dq}{2\pi }}e^{iq\left(f_k\alpha \right)}\left(1{\displaystyle _0^{\mathrm{}}}𝑑fP\left(f\right){\displaystyle \frac{1}{\mathrm{\Omega }}}{\displaystyle d^3r\left\{1\mathrm{exp}\left[iq\frac{U^2\left(r\right)}{f}\right]\right\}}\right)^{N1},`$ (18)
where in the last step we have simply changed variables from $`𝐫_1`$ to $`𝐫=𝐫_1𝐫_k`$. Using the averaging methods detailed in Eqs. (8–10) and the Appendix of Ref. , we see that

$$P\left(f\right)=_{\mathrm{}}^{\mathrm{}}\frac{dq}{2\pi }e^{iq\left(f\alpha \right)}\mathrm{exp}\left[u\sqrt{iq}Q\left(\alpha \right)\right],$$
(19)
where

$$Q\left(\alpha \right)=_0^{\mathrm{}}\frac{df}{\sqrt{f}}P\left(f\right),$$
(20)
and

$$u=\frac{16\pi ^{3/2}}{9\sqrt{3}}\mu ^2\frac{N}{\mathrm{\Omega }}.$$
(21)
Using

$$_{\mathrm{}}^{\mathrm{}}\frac{dq}{2\pi }e^{iqA\sqrt{iq}B}=\frac{1}{2\sqrt{\pi }}\frac{B}{A^{3/2}}\mathrm{exp}\left(\frac{B^2}{4A}\right)\theta \left(A\right),$$
(22)
for $`A`$ and $`B`$ real, we have

$$P\left(f\right)=\frac{1}{2\sqrt{\pi }}\frac{uQ\left(\alpha \right)}{\left(f\alpha \right)^{3/2}}\mathrm{exp}\left[\frac{u^2Q^2\left(\alpha \right)}{4\left(f\alpha \right)}\right]\theta \left(f\alpha \right).$$
(23)
and so

$`Q\left(\alpha \right)`$ $`=`$ $`{\displaystyle _\alpha ^{\mathrm{}}}{\displaystyle \frac{df}{\sqrt{f}}}{\displaystyle \frac{1}{2\sqrt{\pi }}}{\displaystyle \frac{uQ\left(\alpha \right)}{\left(f\alpha \right)^{3/2}}}\mathrm{exp}\left[{\displaystyle \frac{u^2Q^2\left(\alpha \right)}{4\left(f\alpha \right)}}\right]`$ (24)
$`=`$ $`{\displaystyle \frac{1}{\sqrt{\alpha }}}\mathrm{exp}\left[{\displaystyle \frac{u^2Q^2\left(\alpha \right)}{4\alpha }}\right]\mathrm{erfc}\left[{\displaystyle \frac{uQ\left(\alpha \right)}{2\sqrt{\alpha }}}\right].`$ (25)
We also have

$`{\displaystyle \frac{1}{f}}`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{df}{f}}P\left(f\right)`$ (26)
$`=`$ $`{\displaystyle \frac{1}{\alpha }}{\displaystyle \frac{\sqrt{\pi }}{2}}{\displaystyle \frac{uQ\left(\alpha \right)}{\alpha ^{3/2}}}\mathrm{exp}\left[{\displaystyle \frac{u^2Q^2\left(\alpha \right)}{4\alpha }}\right]\mathrm{erfc}\left[{\displaystyle \frac{uQ\left(\alpha \right)}{2\sqrt{\alpha }}}\right],`$ (27)
or

$$\frac{1}{f}=\frac{1}{\alpha }\left[1\frac{\sqrt{\pi }}{2}uQ^2\left(\alpha \right)\right].$$
(28)
Thus we can solve numerically the transcendental relation for $`Q\left(\alpha \right)`$ given in Eq. (25) and plug it into Eq. (28) to find $`1/f`$. It is important to note that Eq. (25) defines $`Q`$ as an analytic function of $`\alpha `$ in the half-plane where $`\mathrm{Re}\alpha >0`$. Thus $`1/f`$ can be analytically continued, in the variable $`\omega =i\alpha `$, from the positive imaginary axis to the entire upper half-plane.

### C $`U`$ and $`V`$ Processes in the Absence of $`\mathrm{\Delta }`$

Now we allow the $`V`$ process to be present, but we still ignore the detuning, $`\mathrm{\Delta }`$. We start with the equations

$`f_0`$ $`=`$ $`\alpha +{\displaystyle \underset{k=1}{\overset{N}{}}}{\displaystyle \frac{V_{0k}^2}{f_k}},`$ (29a)
$`f_k`$ $`=`$ $`\alpha +{\displaystyle \underset{l=1,lk}{\overset{N}{}}}{\displaystyle \frac{U_{kl}^2}{f_l}}.`$ (29b)
Repeating the analysis we did before, we have

$$P_0\left(f_0\right)=\delta \left(f_0\alpha \underset{k=1}{\overset{N}{}}\frac{V_{0k}^2}{f_k}\right)\underset{l=1}{\overset{N}{}}P\left(f_l\right)df_l,$$
(30)
which leads to the form

$$P_0\left(f_0\right)=_{\mathrm{}}^{\mathrm{}}\frac{dq}{2\pi }e^{iq\left(f_0\alpha \right)}\left\{_0^{\mathrm{}}𝑑fP\left(f\right)\frac{1}{\mathrm{\Omega }}d^3r\mathrm{exp}\left[iq\frac{V^2\left(r\right)}{f}\right]\right\}^N.$$
(31)
At this point it is clear from our previous work that this equation leads to

$$P_0\left(f_0\right)=\frac{1}{2\sqrt{\pi }}\frac{vQ\left(\alpha \right)}{\left(f_0\alpha \right)^{3/2}}\mathrm{exp}\left[\frac{v^2Q^2\left(\alpha \right)}{4\left(f_0\alpha \right)}\right]\theta \left(f_0\alpha \right),$$
(32)
where

$$v=\frac{16\pi ^{3/2}}{9\sqrt{3}}\mu \mu ^{}\frac{N}{\mathrm{\Omega }},$$
(33)
and $`Q\left(\alpha \right)`$ is still given by Eq. (25), and depends on $`u`$. Therefore we obtain the result

$`{\displaystyle \frac{1}{f_0}}`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{df_0}{f_0}}P_0\left(f_0\right)`$ (34)
$`=`$ $`{\displaystyle \frac{1}{\alpha }}{\displaystyle \frac{\sqrt{\pi }}{2}}{\displaystyle \frac{vQ\left(\alpha \right)}{\alpha ^{3/2}}}\mathrm{exp}\left[{\displaystyle \frac{v^2Q^2\left(\alpha \right)}{4\alpha }}\right]\mathrm{erfc}\left[{\displaystyle \frac{vQ\left(\alpha \right)}{2\sqrt{\alpha }}}\right],`$ (35)
which again can be analytically continued.

### D $`U`$ and $`V`$ Processes in the Presence of $`\mathrm{\Delta }`$

If we further consider the effect of a detuning, $`\mathrm{\Delta }`$, 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=\mathrm{Re}f_0`$ and $`f_0^I=\mathrm{Im}f_0`$. We have

$$P_0(f_0^R,f_0^I)=\delta \left(f_0^R\alpha \underset{k=1}{\overset{N}{}}\frac{V_{0k}^2}{f_k}\right)\delta \left(f_0^I\mathrm{\Delta }\right)\underset{l=1}{\overset{N}{}}P\left(f_l\right)df_l.$$
(36)
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{\pi }}\frac{vQ\left(\alpha \right)}{\left(f_0^R\alpha \right)^{3/2}}\mathrm{exp}\left[\frac{v^2Q^2\left(\alpha \right)}{4\left(f_0^R\alpha \right)}\right]\theta \left(f_0^R\alpha \right)\delta \left(f_0^I\mathrm{\Delta }\right).$$
(37)
It follows, then, that

$`{\displaystyle \frac{1}{f_0}}`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑f_0^R{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑f_0^I{\displaystyle \frac{1}{f_0^R+if_0^I}}P_0(f_0^R,f_0^I)`$ (38)
$`=`$ $`{\displaystyle \frac{1}{\alpha +i\mathrm{\Delta }}}{\displaystyle \frac{\sqrt{\pi }}{2}}{\displaystyle \frac{vQ\left(\alpha \right)}{\left(\alpha +i\mathrm{\Delta }\right)^{3/2}}}\mathrm{exp}\left[{\displaystyle \frac{v^2Q^2\left(\alpha \right)}{4\left(\alpha +i\mathrm{\Delta }\right)}}\right]\mathrm{erfc}\left[{\displaystyle \frac{vQ\left(\alpha \right)}{2\sqrt{\alpha +i\mathrm{\Delta }}}}\right].`$ (39)
### 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=\alpha +\underset{l=1,}{\overset{N}{}}\frac{U_{kl}^2}{f_k},$$
(40)
which gives explicitly

$$f_k=\frac{1}{2}\left(\alpha +\sqrt{\alpha ^2+4\underset{l}{}U_{kl}^2}\right).$$
(41)
The standard averaging procedure then leads to

$$\frac{1}{f}=\frac{2}{u\sqrt{\pi }}+\frac{\alpha }{u^2}\left[\mathrm{exp}\left(\frac{u^2}{\alpha ^2}\right)\mathrm{erfc}\left(\frac{u}{\alpha }\right)1\right].$$
(42)
Alternatively, we can replace $`1/f_l`$ by its average, obtaining

$$f_k=\alpha +\underset{l=1}{\overset{N}{}}U_{kl}^2\frac{1}{f},$$
(43)
which leads to

$$\alpha \frac{1}{f}+\frac{\sqrt{\pi }}{2}\frac{u}{\sqrt{\alpha }}\sqrt{\frac{1}{f}}\mathrm{exp}\left(\frac{u^2}{4\alpha }\frac{1}{f}\right)\mathrm{erfc}\left(\frac{u}{2\alpha }\sqrt{\frac{1}{f}}\right)=1.$$
(44)
Eqs. (40) and (43) 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. (40) and (43) 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. (33). We have then results as a function of time or frequency and of two parameters: $`\mathrm{\Delta }/v`$, where $`\mathrm{\Delta }`$ is the detuning from resonance, and $`u/v`$, where $`u`$ is the equivalent of $`v`$ for the $`spps`$ transition and is defined in Eq. (21). 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. . 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 ($`NN^{}`$). 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 $`\mathrm{\Delta }=0`$ finding $`a_0`$ is equivalent to finding the on-site Green’s function $`G_{00}`$ for a gas of $`N+1`$ atoms. Then $`\left(1/\pi \right)\mathrm{Im}G_{00}\left(\omega \right)=\left(1/\pi \right)\mathrm{Re}a_0\left(\omega \right)`$ 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. 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\left|a_0\right|^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 $`\mathrm{\Delta }`$ 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. . 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 $`\left(\sqrt{\pi }/2\right)vt`$, while each realization is quadratic in $`t`$, as already discussed in Ref. . It can be seen from the graphs that the linear rise law, $`\left(\sqrt{\pi }/2\right)vt`$, is independent both of $`\mathrm{\Delta }`$ and of $`u`$ and holds for a substantial range of $`vt`$.

### B Results for $`a_0\left(t\right)`$

The four parts of Figs. 2 and 3 display the numerical simulations of $`\mathrm{Re}a_0`$ and $`\mathrm{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 $`\mathrm{\Delta }`$ 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.  for computing $`1\left|a_0\right|^2`$. We give here only the following result for the Laplace transform:

$$a_0\left(i\alpha \right)=\frac{1}{\alpha +i\mathrm{\Delta }}\frac{\sqrt{\pi }}{2}\frac{v}{\left(\alpha +i\mathrm{\Delta }\right)^{3/2}\sqrt{\alpha }}\mathrm{exp}\left[\frac{v^2}{4\alpha \left(\alpha +i\mathrm{\Delta }\right)}\right]\mathrm{erfc}\left[\frac{v}{2\sqrt{\alpha \left(\alpha +i\mathrm{\Delta }\right)}}\right].$$
(45)
From this, $`a_0\left(t\right)`$ has been obtained by expanding in powers of $`1/\alpha `$ and inverting the Laplace transform term by term. As in Ref. , 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_01\frac{\sqrt{\pi }}{2}vti\mathrm{\Delta }t$$
(46)
Looking at the graphs, it seems that the initial decrease of $`\mathrm{Re}a_0`$ depends on $`u`$, but this perception only indicates that the range of validity of the linear approximation (46) is small. On the other hand, the graphs of $`\left|a_0\right|^2`$ in Fig. 1 look linear and independent of $`u`$ over a wide range of $`vt`$. Note also that $`\left|a_0\right|^21\sqrt{\pi }vt`$ decreases initially twice as fast as $`\left|a_0\right|^21(\sqrt{\pi }/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 $`\left|a_0\right|^2`$ is not even a fair approximation to $`\left|a_0\right|^2`$.

### C Results for the Spectral Density

Figs. 4 to 7 present graphs of $`\mathrm{Re}a_0\left(\omega \right)`$ that correspond to the graphs of $`a_0\left(t\right)`$ in Figs. 2 and 3. For $`u=v`$ and $`\mathrm{\Delta }=0`$, $`\mathrm{Re}a_0\left(\omega \right)`$ is the exciton spectral density (times $`\pi `$); 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. (45), obtained from Eqs. (35) and (25) by analytic continuation. We compare spectral densities, rather than time graphs, in part because obtaining $`a_0\left(t\right)`$ 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. (7) as

$$\mathrm{Re}a_0\left(\omega +i\epsilon \right)=\mathrm{Im}\underset{m=0}{\overset{N}{}}\frac{\beta _{0m}\beta _{m0}}{\omega \lambda _m+i\epsilon }$$
(47)
in the limits $`\epsilon 0`$ and $`N\mathrm{}`$. In practice it is difficult to go beyond $`N=1000`$, and one must choose $`\epsilon `$ 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 $`\omega =0`$, and if $`\epsilon `$ 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 $`\epsilon /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 $`\epsilon /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 $`\epsilon `$ 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. Thus, $`\epsilon `$ is almost equivalent to the $`\gamma `$ of the Lorentzian broadening model of Ref. , the difference being that in Ref.  $`i\gamma `$ is added to $`\omega `$ for the Fourier transform of Eq. (1b) only, while here $`i\epsilon `$ is added to every $`\omega `$. As $`\epsilon `$ 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 $`\epsilon `$.

For $`u=0`$ (Fig. 4) the spectral density is quite different from a Lorentzian. At $`\mathrm{\Delta }=0`$ there is a quasi-gap around $`\omega =0`$. As $`\mathrm{\Delta }`$ increases the gap extends from $`\omega =0`$ to $`\omega =\mathrm{\Delta }`$.

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 $`\mathrm{\Delta }=0`$ has the appearance of a single, asymmetric line. At $`\mathrm{\Delta }=v`$ and $`\mathrm{\Delta }=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 $`\mathrm{\Delta }`$, 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 $`\epsilon /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.  is replaced by

$$i\dot{c_k}=V_ka_0i\gamma c_k.$$
(48)
Solving the equations in this limit, one finds a pole at

$$\omega =i\frac{\gamma _kV_k^2}{\gamma ^2+\mathrm{\Delta }^2}+\mathrm{\Delta },$$
(49)
which gives a decreasing width as $`u`$ increases, since $`\gamma `$ is proportional to $`u`$. One cannot simply replace $`_kV_k^2`$ by its average, which does not even exist. Instead, one can Fourier transform the average to find the time dependence

$$\mathrm{exp}\left(i\mathrm{\Delta }t\sqrt{\gamma _{eq}t}\right),$$
(50)
where

$$\gamma _{eq}=\frac{v^2\gamma }{\gamma ^2+\mathrm{\Delta }^2},$$
(51)
and $`v`$ has been defined in Eq. (33). A similar argument was applied in Ref.  to find the time dependence of $`\left|a_0\right|^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 ($`\mathrm{\Delta }=0`$), while the correct spectral density is slightly asymmetric, except for $`u=0`$. Correspondingly, $`a_0\left(t\right)`$ is not purely real, as can be seen in Fig. 3. For $`\mathrm{\Delta }0`$, the spectral density at $`\omega <\mathrm{\Delta }`$ is underestimated by the Cayley approximation.

Clearly, the asymmetry of the spectral density at $`\mathrm{\Delta }=0`$ comes from the terms with $`lm`$ in the denominator of Eq. (9), 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 $`\mathrm{\Delta }=0`$. At finite $`\mathrm{\Delta }`$, the underestimate of $`\mathrm{Re}a_0\left(\omega \right)`$ for $`\omega <\mathrm{\Delta }`$ 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.

### 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. . 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. , where on the basis of a small $`\mathrm{\Delta }`$ 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/v3.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 $`\mathrm{\Delta }=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 $`ppss^{}`$ discussed in Refs.  and  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\left|a_0\right|^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 $`\mathrm{\Delta }`$ shows persistent oscillations with period $`\mathrm{\Delta }`$, in particular when exciton propagation is a strong effect, as for $`u=4v`$. Correspondingly, the spectral density for exciton creation is strongly peaked at $`\mathrm{\Delta }`$, 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 ($`uv`$).

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 $`\left|a_0\right|^2`$, which is notoriously more difficult to obtain than $`a_0`$ and quite different from $`\left|a_0\right|^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. 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\left(r\right)1/r^3`$ is long range and has the distribution $`P\left(U\right)1/U^2`$ for a random distribution of $`r`$. The distribution of the cumulative variable $`𝒰`$, where $`𝒰^2=_kU^2\left(𝐫_k\right)`$, plays a major role in the theory. It is given by

$$P\left(𝒰\right)=\frac{1}{\sqrt{\pi }}\frac{u}{𝒰^2}\mathrm{exp}\left(\frac{u^2}{4𝒰^2}\right),$$
(52)
where $`u`$ is given by the appropriate equivalent of Eq. (21).The angular dependence of $`U`$ modifies the value of $`u`$, leaving unchanged the functional form of Eq. (52). See Ref. , Appendix A. Note that $`P\left(𝒰\right)`$ is not Gaussian and the central limit theorem does not apply, because $`P\left(U\right)`$ 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. (52) makes the Cayley tree equations for $`a_0`$ exactly soluble in this case.

###### Acknowledgements.

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.


text:
# Turbulence in a free surface

## Abstract

We report an experimental and numerical study of turbulent fluid motion in a free surface. The flow is realized experimentally on the surface of a tank filled with water stirred by a vertically oscillating grid positioned well below the surface. The effect of surface waves appears to be negligible so that the flow can numerically be realized with a flat surface and stress-free boundary conditions. The surface flow is unconventional in that it is not incompressible and neither energy nor enstrophy are conserved. Nevertheless, according to both experiment and numerical simulation, the second order structure function $`S_2(R)`$ scales essentially as for a three-dimensional system. However, the surface flow seems to be more intermittent.
PACS numbers: 47.27.Gs, 68.10.-m, 06.30.Gv

Statistically stationary turbulence is readily produced in a tank of water by vertically oscillating a grid below the water surface. The properties of this three-dimensional volume turbulence have been studied extensively . Particles floating on the free surface are driven by the horizontal component of the turbulent flow below. While this is a classic problem, generally studied in an oceanographic context where the “particles” are micro organisms, such as phytoplankton , or Lagrangian drifters , controlled laboratory measurements of the relative velocity of floating particles do not appear to exist, despite a very large number of oscillating grid experiments which probe the fluid flow in the bulk and even near the surface . Such surface flows also seem to have received little attention by those interested in the fundamental aspects of two or three dimensional turbulence . Our experimental and theoretical results demonstrate that these two-dimensional surface flows have unusual properties that make them interesting in connection with scalar transport. The work opens the possibility of studying certain aspects of passive scalar dynamics in divergent flows beyond the ones used for the Kraichnan model which are Gaussian and white in time . Their properties are also needed for a hydrodynamic modelling of the distribution of flotsam driven by underlying turbulence or chaos . Our study is also related to recent work on quasi two-dimensional turbulence, which is probed by tracking the motion of floating particles . The surface flows being comprised of particles constrained to the surface in both the experimental and numerical work, are qualitatively different from the motion of the driving fluid.

Although the fluid is incompressible, the inability of the floaters to enter the bulk flow assures that they will have a non-zero two-dimensional divergence. The compressibility of the floating particles causes them to clump, an affect seen in the numerical and experimental realizations. The main differences between their motion and that of the underlying fluid comes from the possibility that the floaters can exchange kinetic energy and enstrophy with the water molecules below. Hence, the energy and enstrophy of the floaters still would not be constants of the motion, even if the underlying fluid were invisid and undriven . With the above conservation laws absent, one lacks the traditional dimensional arguments to estimate the scaling forms of the velocity and vorticity structure functions, $`S_2(R)`$ and $`S_\omega (R)`$. The structure functions are the main focus of our analysis and are defined as

$$S_2(R)=[(𝐯(𝐱+𝐑)𝐯(𝐱))𝐑/R]^2$$
(1)
and

$$S_\omega (R)=[\omega (𝐱+𝐑)\omega (𝐱)]^2,$$
(2)
where $`𝐑`$ is the separation between two points and both $`𝐱=\{x,y\}`$ and $`𝐑`$ are in a horizontal plane. The vorticity, of magnitude $`\omega (x,y)`$, is perpendicular to that plane. A third quantity characterizing the surface flow is the dimensionless compressiblity coefficient $`𝒞`$, defined as

$$𝒞=(𝐯)^2/(𝐯)^2,$$
(3)
which lies between 0 and 1 if the turbulence is isotropic . As usual, $`\mathrm{}`$ denotes an ensemble average.

The measurements were made in a square plexiglass tank filled with water to a height that was varied. The bars of the grid were PVC and were square in cross section with dimensions of 1 cm. The grid is also square and conforms closely with the square shape of the plexiglass box. To minimize flexing of the grid when it is oscillating, a number of the PVC bars are replaced by brass bars. For some of the measurements, a vertical post was placed at the center of the square grid to further suppress flexing, thereby reducing the amplitude of the surface waves. The apertures of the grid are also square, their dimensions being $`L`$=3.66 cm. The grid is supported by four vertical rods which run through pillow blocks, assuring its smooth vertical motion. All the parameters used in our set-up were typical of prior oscillating grid experiments , except that the lateral dimensions of the tank were somewhat larger.The grid is driven sinusoidally from above by a 1/4 hp motor. The vertical amplitude of grid motion was $`A=`$ 1.25 cm. Whereas these parameters were varied by only a small amount, the distance $`Z`$ between the grid and the water surface (at $`z_0`$=0) was rather widely varied. Reducing $`Z`$ increases the Reynolds number of the turbulence on the surface and in the bulk close to the surface ($`z_0>0`$). For all the measurements presented here, the oscillation frequency $`f`$ of the grid was 4.5 Hz $`\pm `$ 0.3 Hz.

The flow was measured using the technique of particle imaging velocimetry (PIV). There one illuminates the seed particles with a (horizontal) sheet of laser light and tracks their motion photographically. A CCD camera was centered at a distance that was 23 cm from the nearest wall. It was mounted above the tank where we could capture images of the surface particles or those in the interior. The particle illumination was furnished by a Nd-Yag laser, with the beam spread into a horizontal sheet by a cylindrical lens. This sheet of light was typically 10 cm wide and 0.5 cm thick. The floaters were particles of various diameters ranging from 10 to 50 $`\mu `$m. The interior measurements were made using neutrally buoyant polystyrene spheres of diameter 10 $`\mu `$m. The images were acquired using a commercial particle imaging velocimetry apparatus  and the particle tracking was done using in-house software. To form a statistically reliable ensemble average, each measurement of $`S_2(R)`$ lasted several minutes and included $``$ 200 image pairs. Each $`S_2(R)`$ was determined from roughly a million particle pairs.

FIG. 1. Log-log plot of the second moment of the (longitudinal) velocity difference at the indicated distances ($`Z`$) between the surface and the oscillating grid. The oscillating frequency of the grid and the surface-to-grid distances, $`Z`$, are indicated. Other parameters are $`Z=8.2`$ cm: $`Re_\lambda `$=80, $`l_0`$=3.6 cm, $`v_{rms}`$=1.27 cm/s; $`Z`$= 18.2 cm: $`Re_\lambda `$=50, $`l_0`$=1.5 cm, $`v_{rms}`$=0.7 cm/s. For both sets of measurements $`𝒞`$=0.5 $`\pm `$ 0.1.

The surface and bulk turbulence are characterized by the outer scale of the turbulence $`l_0`$ and the Taylor microscale $`\lambda `$. The outer scale is defined as the integral of the normalized velocity autocorrelation function $`C(R)=v(x+R)v(x)/v_{rms}^2`$, where the brackets indicate an average over points $`x,y`$ and over many images. Here $`v_{rms}^2(v_i(x,y)v_i(x,y))^2`$ where $`i`$ is one of the two arbitrarily chosen directions in the horizontal plane. Separate experiments established that the turbulence is almost isotropic in the horizontal plane (see below). The Taylor microscale is defined here as $`\lambda =\sqrt{(v_{rms}^2)/(dv(x)/dx)^2}`$. The Taylor microscale Reynolds number is $`Re_\lambda =v_{rms}\lambda /\nu `$, where $`\nu `$ is the kinematic viscosity (in water $`\nu `$=0.01 cm2/s). The surface motion is further characterized by the vertical component of the vorticity $`\omega `$ and by the dimensionless compressibility coefficient $`𝒞`$ \[see Eq. (3)\].

Figure 1 is a log-log plot of $`S_2(R)`$ in the surface measured for $`Z`$= 18.2 cm (closed circles) and for $`Z`$= 8.2 cm (open circles). The grid oscillation frequency for the two sets of measurements are indicated. The outer scale of the turbulence $`l_0`$ at $`Z`$=18.2 cm and 8.2 cm are 1.5 and 3.6 cm, respectively. The Taylor microscale Reynolds numbers $`Re_\lambda `$ are 80 ($`Z`$ = 8.2 cm) and 50 ($`Z`$= 18.2 cm). Technical limitations blocked measurement of $`S_2`$ at $`R`$ less than 1 mm. The measurements in Fig. 1 suggest that $`S_2(R)`$ is approaching its saturation value $`2v_{rms}^2`$ at large $`R`$. The more interesting case is $`S_2`$ at smaller $`R`$, where one finds that $`S_2(R)R^\zeta `$, with $`\zeta `$=1.6 and 1.4 at $`Z`$= 8.2 and 18.2 cm, respectively. The dissipative range at even smaller $`R`$, where $`S_2(R)R^2`$, was not accessible in the laboratory experiment . For these, and the other surface measurements, the dimensionless compressibility $`𝒞`$ defined above, was 0.5 $`\pm `$ 0.1. Reynolds numbers in the experiment and the simulation are too small to develop a significant interval with algebraic scaling in $`S_2(R)`$ between the dissipative, $`S_2(R)R^2`$, and the saturation range $`Rl_0`$, where $`S_2(R)R^0`$. Therefore the extended self similarity (ESS) representation is used  and a local scaling exponent $`D_{2,\mathrm{\hspace{0.17em}3}}(R)=d\mathrm{log}_eS_2(R)/d\mathrm{log}_eG_3(R)`$ for all $`R`$ in the interval 0.1 cm $`<R<`$ 7 cm is found to give $`0.70\pm 0.10`$ (see upper panel of Fig. 2). For the ESS plots, the third order structure function $`G_3(R)=\left|(𝐯(𝐱+𝐑)𝐯(𝐱))𝐑/R\right|^3`$ is needed.

Since the surface is entirely driven by the grid-generated turbulence in the water, it is important to measure $`S_2(R)`$ below the surface. We have made such measurements at distances $`z_0`$ in the range 0.5 to 2 cm and for various grid-to-surface distances $`Z`$. At $`Z`$ 9.6 cm, $`z_0`$= 2.0 cm, the second moment is very well fitted to a power law form in the interval 0.1 cm $`<R<`$ 4 cm. In this interval $`\zeta `$=0.9 $`\pm `$ 0.2, which is somewhat larger than the Kolmogorov value of 2/3 . Other relevant parameters for this bulk measurement are $`Re_\lambda `$= 100 and $`l_0`$= 3.1 cm. The corresponding parameters measured on the surface are $`Re_\lambda `$= 40 and $`l_0`$= 3.5 cm. At all values of $`R`$, the bulk value of the second moment is roughly an order of magnitude larger than $`S_2(R)`$ measured at the surface. A dimensional argument suggests the crossover between bulk and surface behavior should take place at a depth $`\delta z`$ below the surface given by $`\delta z\nu /\stackrel{~}{v}`$ where $`\stackrel{~}{v}`$ is the bulk rms value of $`v_z`$. This gives $`\delta z`$ 0.1 mm. The numerical studies  indicate that the prefactor in this dimensional estimate can be as large as 10, making $`\delta z`$ 1 mm, which is still too close to the surface for us to measure.

FIG. 2. Local scaling exponent, $`D_{2,3}(R)`$, of an extended self similarity plot of second order moment vs. absolute third order moment. Upper panel: The parameters for this experiment are the same as in Fig. 1 at $`Z`$=18.2 cm. Lower panel: Comparison of $`D_{2,3}(R)`$ in the bulk at $`z_0`$=0.5 (triangles) with that on the two surfaces at $`z_0`$=0 and 1 (diamonds). Data were taken from a numerical simulation at $`Re_\lambda `$ 100. Dashed lines in both panels are drawn at $`D_{2,3}(R)`$=0.73 and 0.71, respectively.

We consider the possibility that the surface velocity fluctuations originate from gravity and capillary waves generated by the 4.4 Hz oscillation frequency of the grid. Even though the amplitude of the surface waves was rather small ($`A`$ 0.5 mm) it is comparable to the thickness of the boundary layer near the surface. The appropriate wave turbulence theory treats the nonlinear coupling of surface modes $`\mathrm{\Omega }(k)`$ of different wavenumbers $`k`$ . This theory is of Hamiltonian form, which requires that the (2D) velocity field in the theory is obtained from a velocity potential, so that the vorticity of the surface particles is zero. Measurements of the (vertically-directed) vorticity of the surface particles show that it is quite large. Figure 4 is a log-log plot of the second moment of the surface vorticity fluctuations, $`S_\omega (R)`$, for a grid-to-surface distance $`Z`$= 8.5 cm. This function saturates at a rather small value of $`R`$ and approaches the limiting value $`2\omega ^2`$. A dimensionless measure ($`Q`$) of the surface vorticity is $`S_\omega (R)R^2/S_2(R)`$. At $`R`$ equal to the outer scale of the turbulence, $`l_03`$ cm, $`Q`$ 1, which suggests that it is probably too large to satisfy the zero-vorticity requirement of wave turbulence theory, even though this theory successfully accounts for the motion of particles floating on a container of vertically oscillating fluid where capillary waves are excited .

The experimental observations are, however, in good agreement with a model for turbulence in a flat surface bounding a three-dimensional volume with fully developed turbulence governed by the incompressible Navier-Stokes equation . This fluid is in a computational box with periodic boundary conditions applied to the lateral edges but with the fluid particles on the upper and lower surfaces obeying free–slip boundary conditions, with $`v_z=_zv_x=_zv_y=0`$. The equations were solved by means of a pseudospectral method with resolution $`N_x\times N_y\times N_z`$ of $`256\times 256\times 65`$. A statistically steady state of the turbulent flow was achieved by forcing in a wavenumber range comparable to the inverse of the box height (for more details see ). Here $`Re_\lambda `$ 100.

FIG. 3. Deviations of the local extended self-similarity (ESS) slope from the classical Kolmogorov scaling for different heights $`z_0`$ of the averaging plane and for different orders $`n=2`$ to $`n=6`$. All data points are calculated from $`2\times 114`$ statistically independent samples at $`Re_\lambda `$ 100 and by averaging $`D_{n,\mathrm{\hspace{0.17em}3}}(R,z_0)`$ for scales $`R`$ between $`18\eta `$ and $`41\eta `$.

The local ESS scaling exponent $`D_{2,\mathrm{\hspace{0.17em}3}}(R)`$ lies in the range between 0.67 and 0.74 in the interval 5 $`<R/\eta <`$ 100 with the Kolmogorov length $`\eta =(\nu ^3/ϵ)^{1/4}`$ (lower panel of Fig. 2) . This result is in good agreement with the experimental data also shown in Fig. 2 (upper panel). The simulation results show a plateau between $`18\eta `$ and $`41\eta `$ for the local scaling exponents in both cases, at the surface and in the bulk. Intermittency corrections in the surface seem to be stronger than in the bulk. In Fig. 3, a systematic investigation of these deviations from the classical Kolmogorov scaling with respect to $`z_0`$, the vertical position of the horizontal plane of averaging is presented. Our results indicate increasing deviations toward the free surface. The two-dimensional surface compressibility $`𝒞`$ in the calculation is 0.5, in excellent agreement to the laboratory-measured value of $`0.5\pm 0.1`$. This value of $`𝒞`$ is expected for isotropic turbulence.

FIG. 4. Second moment of vorticity differences for the floating particles, at $`Z`$=9.6 cm. Here $`f`$=4.7 Hz, $`Re_\lambda `$=35, $`l_0`$=1.7 cm, and $`v_{rms}`$=0.67 cm/s. Again $`𝒞`$ 0.5. Inset: Model calculation of $`S_\omega (R)`$ at the free-slip surfaces, with $`Re_\lambda `$ 100. Here $`𝒞`$=0.5.

The above model was also applied to the calculation of $`S_\omega (R)`$, and the results are presented in the inset of Fig. 4. Both the calculated and measured $`S_\omega (R)`$ are qualitatively similar and show that this function rises rapidly with increasing $`R`$ in these log-log plots and reaches its saturation value, 2$`\omega ^2`$ at a rather small value of $`R`$ .

In summary, we have described a laboratory experiment and a model calculation for particles that float on a sea of three-dimensional turbulence. Even though the surface particles move in a plane, there appears to be no basis for describing the floating particle system in terms of the laws of two-dimensional turbulence or of turbulent surface waves. Important characteristics of the floating particle system are its large compressibility, large vorticity and the absence of conservation laws for energy and enstrophy. The measurements and calculations presented here indicate the existence of an inertial range where $`S_2(R)`$ scales as $`R^\zeta `$, with $`\zeta `$ quite close to 2/3, as in 3D flows and in the inverse cascade of 2D turbulence.

In this work we have profited from interactions with P. Alstrøm, S. Banerjee, D. Lohse, C. Cheung, G. Falkovich, V. Horváth, M. Rivera, E. Schröder, X. L. Wu, and M. Vergassola. The particle tracking software of M. Rivera was very important to the success of this experiment. Derek Hohman provided valuable assistance with the experiment. The work has been strongly influenced by the attendance of two of us (W.I.G., B.E.) at the Institute for Theoretical Physics at Santa Barbara. Research support from NASA, the National Science Foundation, and the European Community is gratefully acknowledged. We thank the John von Neumann–Institut für Computing at the Forschungszentrum Jülich for support and computing time on a Cray T-90 for the numerical simulations.