text:
# Introduction

## Introduction

Parton cascade codes that take a space-time approach to model the microscopic processes 111Prominent examples are VNI and pcpc; cf. the OSCAR archive at
http://nt3.phys.columbia.edu/people/molnard/OSCAR/ .  of an ultrarelativistic heavy-ion interaction are – in spite of their QCD bells and whistles – by necessity based on some kind of *classical many-particle dynamics*.

### Problems:

* space-time-based cascade models suffer from the consequences of the

No-Interaction Theorem (NIT) \[CJS63\]:

The only consistent many-particle Hamiltonian theory that is Poincaré-covariant is that of a system of free particles.
* the procedure to determine the sequence of the binary parton interactions (“SBPI”) – an essential aspect of space-time-based cascade models – is by necessity an artificial and ad-hoc feature of these codes.

The only way to circumvent the NIT is to loosen its assumptions:

1. forget about Poincaré covariance $``$ VNI
In this approach the SBPI depends on the initially chosen frame of reference. Einstein causality remains a problem.
2. introduce a many-times formalism, e.g. by formulating the model in $`8N`$-dimensional phase space ($`N`$ is the particle number) 222We have shown previously that for $`N=2`$ all known Poincaré-covariant formalisms are equivalent to ours (cf. \[PNB94\]).
$``$ pcpc
The Poincaré covariance of this model seems to guarantee Einstein causality; the SBPI of the code, however, deserves closer scrutiny.

## Basic Covariant Structure of pcpc

pcpc is a hybrid of a *classical* dynamics approach that governs the evolution of the system between binary parton interactions, and a parton interaction model with QCD ingredients (parton distribution functions, pQCD cross sections and ‘DGLAP evolution’) 333For details, cf. \[PNB94, BMGMN00\].

* The dynamical evolution of the system is parametrized by a Poincaré scalar, $`s`$ (in contrast to the usual $`t`$, the proper time of an external observer, as measured in some external frame)
* The phase-space variables of the $`N`$ partons (quarks, gluons) are covariant 4-vectors $`x_{i}^{}{}_{}{}^{\mu }(s)`$, $`p_{i}^{}{}_{}{}^{\mu }(s)`$ . Between binary interactions $`N`$ is fixed, but can (and does) vary due to parton creation in the QCD-governed parton interactions
* The ($`8N`$-dimensional) interaction term of the Poincaré-*invariant* ‘Hamiltonian’ depends on the following Poincaré-invariant ‘4-distances’

$$d_{ij}^{}{}_{}{}^{2}(s):=\widehat{x}^2=(x\widehat{x})=(\widehat{x}x)$$

$$\widehat{x}^\mu :=x^\mu \frac{(xp)}{p^2}p^\mu $$

\[ $`x^\mu (s):=x_{i}^{}{}_{}{}^{\mu }(s)x_{j}^{}{}_{}{}^{\mu }(s),p^\mu (s):=p_{i}^{}{}_{}{}^{\mu }(s)+p_{j}^{}{}_{}{}^{\mu }(s)`$ are the relative 4-distance and total 4-momentum of particles $`i`$ and $`j`$ \]

The physical reasons why in Poincaré-covariant dynamics the interactions can only depend on these “orthogonal projections” $`\widehat{x}^\mu `$ have been given before \[PNB94\].
* Binary parton interactions occur at the $`s`$ determined by $`d_{ij}^{}{}_{}{}^{2}(s)`$ being at a minimum 444The interaction term in the Hamiltonian thus can be thought of as a sum of $`\delta `$-functions.  (in the cms of partons i, j, $`d_{ij}^{}{}_{}{}^{2}`$ is the minimal 3-distance of approach)
* Between interactions, partons move along free trajectories:

$$x_{i}^{}{}_{}{}^{\mu }(s)=\frac{p_{i}^{}{}_{}{}^{\mu }(s)}{m_i}(ss_0)+x_{i}^{}{}_{}{}^{\mu }(s_0)$$

## Logical flow in pcpc

In contrast to most other cascade codes, the ‘time step’ in pcpc is not a computational artifice, but is given by the formalism itself:

* if there were only 2 partons, their interaction would occur at the value of the evolution parameter $`s`$ at which they would reach their minimal $`d_{ij}^{}{}_{}{}^{2}(s)`$ if they were indeed alone in the world,
* so a ‘time table’ is kept, containing, for every pair $`i,j`$ of partons, its potential minimal approach $`d_{ij}^{}{}_{}{}^{2}(s)`$ (i.e. their minimal approach if $`i,j`$ were the only partons in the whole system), and the corresponding $`s`$,
* then this table is searched for the *smallest* (‘earliest’) $`s`$. At this $`s`$ the potential interaction of the corresponding pair *will actually occur* (because all other potential interactions are ‘later’),
* the interaction of the pair $`i,j`$ will change the world lines of these two partons (and possibly create further partons). Therefore, the time table is updated, with new values of $`d_{ij}^{}{}_{}{}^{2}(s)`$ for all pairs involving either parton $`i`$ or parton $`j`$ (or the newly created partons); and the code loops.

Thus the code follows exactly the sequence of binary parton interactions as parametrized by the monotonically increasing *Poincaré-invariant* evolution parameter $`s`$. The logic is summarized in the following flow diagram:

## Do the parton interactions preserve Einstein causality?

Prima facie this does not seem to be so: the 4-vector $`x_{i}^{}{}_{}{}^{\mu }x_{j}^{}{}_{}{}^{\mu }`$ is space-like, and so no signal can be transmitted between events $`x_i`$ and $`x_j`$. But this argument is fallacious: while it would be correct in the frame work of a $`6N`$-dimensional phase space formalism and physical observer time, it does not follow in a many-times formalism with an $`8N`$-dimensional phase space.

Furthermore, it must be remembered that the representation of the (QCD) physics of the heavy-ion reaction in terms of cross sections, parton distribution functions etc. is tantamount to the description of intrinsically quantum processes in a *classical* terminology. But as long as we refrain from trying to look *inside* an individual binary parton interaction with classical concepts, Einstein causality is not infringed upon by such a model (for details on this point cf. \[PNB94\] ).

## Is Einstein causality preserved between binary interactions?

In terms of the cascade picture, a heavy-ion reaction can be seen as one or several disjoint graphs of connected particle world lines. In every connected subgraph, the Poincaré-covariant dynamics guarantees Einstein causality.

Separate disjoint subgraphs, however, can have no causal connection. But by construction, no signals are transmitted between them. Imagining a full quantum-mechanically description of such a system, unconnected subgraphs would correspond to subamplitudes which would simply be multiplied in obtaining the total amplitude.

## Are the initial parton *positions* critically important?

There is no physically convincing argument for how to set the time components of the initial parton 4-vectors $`a_i:=x_{i}^{}{}_{}{}^{\mu =0}(s=s_0)`$; so these must be fixed phenomenologically with some arbitrary prescription. In pcpc they are all set to zero.

Does this imply the choice of a particular initial frame of reference, thus invalidating covariance, or spoiling the Einstein causality of the model? The answer is

NO!

To see this, suppose that we *do not* set the initial $`a_i=0`$, but retain them as free (arbitrary) parameters. We would then find the minimum of the $`d_{ij}^{}{}_{}{}^{2}`$ to be formally dependent on the $`a_i`$. But since for any two time-like 4-vectors the invariant quantity $`(x\widehat{y})`$ is simply $`\text{}\stackrel{}{x}\stackrel{}{y}|_{\text{cms}}`$, we find that the $`d_{ij}^{}{}_{}{}^{2}`$ are actually *independent of the* (time components of the) *initial parton positions*. It follows that the parton interactions are implemented in a Poincaré-covariant way, even though their sequence (SBPI) *does* depend on how the initial parameters are chosen.

## Conclusions

To sum up, causality is *not an issue* in discussing the physical validity of a cascade model that uses a space-time approach to the dynamical evolution (provided the model is Poincaré-covariant).

It is, however, important to realize that in such codes the sequence of binary parton interactions (SBPI) is to be considered an essential *part of the model*, and that it is *necessarily phenomenological* in character. In Poincaré-covariant cascade codes, such as pcpc, the SBPI (although remaining a phenomenological prescription) is independent of choice of the frame of reference in which the code is run.

The difficulties with Einstein causality incurred by non-covariance of the SBPI have been discussed many years ago \[KDCDN84\]. In contrast to the situations described in that paper, Einstein causality is preserved in a Poincaré-covariant model such as pcpc, both for the individual binary interactions and for the dynamic evolution of the system as a whole.


format: text

text:
# Jahn-Teller Distortions and the Supershell Effect in Metal Nanowires

## Abstract

A stability analysis of metal nanowires shows that a Jahn-Teller deformation breaking cylindrical symmetry can be energetically favorable, leading to stable nanowires with elliptic cross sections. The sequence of stable cylindrical and elliptical nanowires allows for a consistent interpretation of experimental conductance histograms for alkali metals, including both the electronic shell and supershell structures. It is predicted that for gold, elliptical nanowires are even more likely to form since their eccentricity is smaller than for alkali metals. The existence of certain metastable “superdeformed” nanowires is also predicted.

The effects of shell filling on the abundance spectra of metal clusters have been known for years and explain the existence of clusters with “magic numbers,” corresponding to full electronic shells, which are observed more frequently than others Heer (1996). More recently, electronic shell and supershell structures in alkali metal nanowires have been reported by the Leiden group Yanson et al. (1999, 2000, 2001). The conductance $`G`$ was recorded during the breaking of nanowires in a mechanically controllable break junction (MCBJ), and histograms were built out of thousands of measurements. When the temperature is large enough to allow the wire to explore the phase space of possible shapes, the shell structure shows up as peaks in the histogram of $`\sqrt{G}`$ that are equally spaced, while the supershell structure manifests itself as a modulation of the amplitude of the peaks (See Agraït *et al.* Agraït et al. (2003) for a recent review).

The electron-shell structure in alkali nanowires can be understood qualitatively using the free-electron model of a cylinder Yanson et al. (1999). However, more detailed stability analyses of axisymmetric wires Kassubek et al. (2001); Zhang et al. (2003); Urban and Grabert (2003) have revealed a sequence of stable “magic” radii, with characteristic gaps that are not fully consistent with the nearly perfect periodicity of the experimentally observed peak positions. The deviations can be accounted for neither by the inclusion of disorder Bürki and Stafford (2001), nor by the use of more elaborate self-consistent jellium models Puska et al. (2001). Since gaps in the sequence of cylindrical nanowires arise from a degeneracy of conductance channels, it is natural to assume that a Jahn-Teller deformation breaking the symmetry can lead to more stable deformed configurations.

While the experimental manifestations of electron-shell structure are similar in metal clusters and nanowires, the Jahn-Teller effect plays out quite differently in these two systems. Metal clusters and nanowires differ fundamentally in that surface effects tend to stabilize clusters, while they lead to a Rayleigh instability in nanowires Kassubek et al. (2001). As a result, Jahn-Teller deformations of clusters are very common and typically rather small Bulgac and Lewenkopf (1993); Schmidt et al. (1999) while, as we show here, they only occur for a minority of nanowires and can be rather large. Nevertheless, there are several stable nanowires with elliptic cross-sections (referred to as *elliptical nanowires* henceforth), leading to a sequence of stable cylindrical and elliptical nanowires that allows for a satisfying interpretation of the experiments on electronic shell and supershell structure Yanson et al. (1999, 2000, 2001). In addition, our theory predicts the existence of short-lived superdeformed wires.

We use a nanoscale free-electron model, treating the electrons as a non-interacting Fermi gas confined within the wire by hard-wall boundary conditions Stafford et al. (1997). This continuum model is especially suitable for alkali metals, but can also be applied to other monovalent metals. Macroscopic arguments Zhang et al. (2003) suggest that the surface tension $`\sigma _s`$ should drive a crossover from a crystalline solid to a plastic or fluid below a critical radius $`R_{\mathrm{min}}=\sigma _s/\sigma _Y`$, where $`\sigma _Y`$ is the yield strength. Typically, $`R_{\mathrm{min}}10\text{nm}`$ in metals Zhang et al. (2003). This implies that electronic effects should dominate over atomistic effects for sufficiently small radii. Indeed, a crossover from atomic-shell to electron-shell effects with decreasing radius is observed in both metal clusters Martin (1996) and nanowires Yanson et al. (2001); Díaz et al. (2003); Medina et al. (2003); Mares et al. (2004), justifying a posteriori the use of the nanoscale free-electron model in the later regime.

Guided by evidence of the existence of non-spherical clusters Bulgac and Lewenkopf (1993); Schmidt et al. (1999), we consider straight elliptical nanowires foo (a) aligned along the $`z`$-axis. The wire cross-section is characterized by the ellipse’s two major semi-axes $`a`$ and $`b`$, or equivalently, by the area parameter $`\sigma ^2=ab`$ and the aspect ratio $`\epsilon =a/b`$. When a small perturbation, written as a Fourier series, is added to a wire of length $`L`$, the geometry is characterized by the two functions

$`\sigma (z)`$ $`=`$ $`\overline{\sigma }+\lambda {\displaystyle \underset{q}{}}\sigma _qe^{iqz}\overline{\sigma }+\lambda \delta \sigma (z),`$
$`\epsilon (z)`$ $`=`$ $`\overline{\epsilon }+\lambda {\displaystyle \underset{q}{}}\epsilon _qe^{iqz}\overline{\epsilon }+\lambda \delta \epsilon (z),`$ (1)
where the dimensionless parameter $`\lambda `$ sets the size of the perturbation. The perturbation wave vector $`q`$ must be an integer multiple of $`2\pi /L`$ and, since $`\sigma (z)`$ and $`\epsilon (z)`$ are real, we have $`\sigma _q=\sigma _q^{}`$ and $`\epsilon _q=\epsilon _q^{}`$. The deformation must fulfill another constraint which comes from the fact that, depending on material parameters, the deformed wire tries to find a compromise between a volume conserving deformation and one ensuring electroneutrality Stafford et al. (1999). Here we use the general constraint

$$𝒩k_F^3𝒱\eta (3\pi k_F^2𝒜/8)=\mathrm{const},$$
(2)
where $`𝒱`$ is the volume of the wire, $`𝒜`$ its surface area, and $`k_F=2\pi /\lambda _F`$ the Fermi wavevector. The parameter $`\eta `$ can be adjusted so as to fix the value of the effective surface tension to the material-specific value. In particular, $`\eta =0`$ corresponds to a constant-volume constraint, and $`\eta =1`$ is the constraint of constant (zero temperature) Weyl charge.

A nanowire is an open system, with electrons being injected and absorbed by the leads, so that we need to calculate its grand canonical potential $`\mathrm{\Omega }`$ in order to determine the energetic cost of a deformation. A fully quantum mechanical stability analysis of cylindrical nanowires Urban and Grabert (2003) shows that, if the wires are short enough, and/or the temperature is not too low, their stability is essentially determined by their response to long wavelength perturbations. For these we can use the adiabatic approximation and decouple the transverse and longitudinal motions yielding a series of effectively one-dimensional problems, which can be solved with the WKB approximation, in which $`\mathrm{\Omega }`$ is given by Stafford et al. (1997)

$`\mathrm{\Omega }[T,L;\sigma ,\epsilon ]`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}\text{d}E\left({\displaystyle \frac{f}{E}}\right)\mathrm{\Xi }[E,L;\sigma ,\epsilon ],`$ (3)
$`\mathrm{\Xi }[E,L;\sigma ,\epsilon ]`$ $`=`$ $`{\displaystyle \frac{8E_F}{3\lambda _F}}{\displaystyle _0^L}\text{d}z{\displaystyle \underset{\nu }{}}\left({\displaystyle \frac{EE_\nu (\sigma ,\epsilon )}{E_F}}\right)^{3/2}.`$ (4)
Here $`f=\left(1+\mathrm{exp}[(EE_F)/k_BT]\right)^1`$ is the Fermi function at temperature $`T`$, $`E_F`$ is the Fermi energy, and the $`E_\nu `$’s are the transverse eigenenergies. The sum in Eq. (4) runs over all open channels $`\nu `$, for which $`E_\nu (\sigma ,\epsilon )<E`$. A factor of 2 accounts for spin degeneracy. Note that at zero temperature Eq. (3) reduces to $`\mathrm{\Omega }[T=0,L]=\mathrm{\Xi }[E_F,L]`$.

The transverse eigenenergies can be written as $`E_\nu (\sigma ,\epsilon )=E_F\frac{\gamma _\nu (\epsilon )^2}{(k_F\sigma )^2},`$ where the $`\gamma _\nu `$’s are the zeros of the modified Mathieu functions van den Broek and Peeters (2001), and can be computed numerically. In the limit $`\epsilon 1`$, we recover the result for a cylindrical wire, where most eigenenergies are degenerate. As shown in Fig. 1, this degeneracy is lifted for $`\epsilon 1`$ as a result of the decrease in symmetry, so that an elliptic nanowire can sometimes be more stable than the axisymmetric one, leading to a Jahn-Teller deformation.

The energetic cost of a small deformation of a straight elliptical nanowire can be calculated by expanding Eq. (3) as a series in the parameter $`\lambda `$,

$$\mathrm{\Omega }=\mathrm{\Omega }^{(0)}+\lambda \mathrm{\Omega }^{(1)}+\lambda ^2\mathrm{\Omega }^{(2)}+𝒪(\lambda ^3).$$
(5)
A nanowire $`(\overline{\sigma },\overline{\epsilon })`$ is energetically stable at temperature $`T`$ if $`\mathrm{\Omega }^{(1)}(\overline{\sigma },\overline{\epsilon })=0`$ and $`\mathrm{\Omega }^{(2)}(\overline{\sigma },\overline{\epsilon })>0`$ for every possible deformation $`(\delta \epsilon ,\delta \sigma )`$ satisfying $`𝒩=0`$.

The expansion of Eq. (4) yields

$`{\displaystyle \frac{\lambda _F}{L}}\mathrm{\Xi }^{(0)}`$ $`=`$ $`{\displaystyle \frac{8}{3\sqrt{E_F}}}{\displaystyle \underset{\nu }{}}[E\overline{E}_\nu ]^{3/2},`$ (6)
$`{\displaystyle \frac{\lambda _F}{L}}\mathrm{\Xi }^{(1)}`$ $`=`$ $`{\displaystyle \frac{8}{\sqrt{E_F}}}{\displaystyle \underset{\nu }{}}\overline{E}_\nu \sqrt{E\overline{E}_\nu }\left({\displaystyle \frac{\overline{\gamma }_\nu ^{}}{\overline{\gamma }_\nu }}\epsilon _0{\displaystyle \frac{\sigma _0}{\overline{\sigma }}}\right),`$ (7)
$`{\displaystyle \frac{\lambda _F}{L}}\mathrm{\Xi }^{(2)}`$ $`=`$ $`{\displaystyle \frac{4}{\sqrt{E_F}}}{\displaystyle \underset{q}{}}\left({\displaystyle \genfrac{}{}{0pt}{}{\sigma _q/\overline{\sigma }}{\epsilon _q}}\right)^{}\left(\begin{array}{cc}A_{\sigma \sigma }& A_{\sigma \epsilon }\\ A_{\sigma \epsilon }& A_{\epsilon \epsilon }\end{array}\right)\left({\displaystyle \genfrac{}{}{0pt}{}{\sigma _q/\overline{\sigma }}{\epsilon _q}}\right),`$ (8)
where $`\overline{E}_\nu =E_\nu (\overline{\sigma },\overline{\epsilon })`$, and $`\overline{\gamma }_\nu =\gamma _\nu (\overline{\epsilon })`$. The prime indicates differentiation with respect to $`\epsilon `$. The elements of the matrix $`A`$ in Eq. (8) are given by

$`A_{\sigma \sigma }`$ $`=`$ $`{\displaystyle \underset{\nu }{}}\overline{E}_\nu \left(3\sqrt{E\overline{E}_\nu }{\displaystyle \frac{\overline{E}_\nu }{\sqrt{E\overline{E}_\nu }}}\right),`$
$`A_{\sigma \epsilon }`$ $`=`$ $`{\displaystyle \underset{\nu }{}}\overline{E}_\nu [\{[{\displaystyle \frac{\overline{\gamma }_\nu ^{}}{\overline{\gamma }_\nu }}]^2+{\displaystyle \frac{\overline{\gamma }_\nu ^{\prime \prime }}{\overline{\gamma }_\nu }}\}\sqrt{E\overline{E}_\nu }[{\displaystyle \frac{\overline{\gamma }_\nu ^{}}{\overline{\gamma }_\nu }}]^2{\displaystyle \frac{\overline{E}_\nu }{\sqrt{E\overline{E}_\nu }}}],`$
$`A_{\epsilon \epsilon }`$ $`=`$ $`{\displaystyle \underset{\nu }{}}\overline{E}_\nu {\displaystyle \frac{\overline{\gamma }_\nu ^{}}{\overline{\gamma }_\nu }}\left(2\sqrt{E\overline{E}_\nu }{\displaystyle \frac{\overline{E}_\nu }{\sqrt{E\overline{E}_\nu }}}\right).`$ (9)
The constraint $`𝒩=0`$ on allowed deformations restricts the number of independent Fourier coefficients in Eq. (Jahn-Teller Distortions and the Supershell Effect in Metal Nanowires). Hence we can express $`\sigma _0`$ in terms of the other coefficients and expand it as a series in $`\lambda `$,

$$\sigma _0=\sigma _0^{(0)}(\overline{\epsilon },\overline{\sigma },\epsilon _0)+\lambda \sigma _0^{(1)}(\overline{\epsilon },\overline{\sigma },\{\epsilon _q,\sigma _q\})+𝒪(\lambda ^2),$$
(10)
and eliminate $`\sigma _0`$ from Eqs. (7) and (8), modifying the expressions for $`\mathrm{\Xi }^{(1)}`$ and the stability matrix $`A`$. At zero temperature, the modified condition $`\stackrel{~}{\mathrm{\Xi }}^{(1)}=0`$ determines which wires are stationary states, while the positivity of the modified matrix $`\stackrel{~}{A}`$ reveals the stability of the wire. The results for finite temperature can be derived analogously by evaluating Eq. (3) numerically.

Figure 2 shows the stable geometries as a function of temperature and conductance, with the surface tension adjusted to the value $`0.22`$ N/m (i.e. $`\eta =0.93`$), corresponding to Na Tyson and Miller (1977). The color represents the value of the aspect ratio $`\epsilon `$ and the $`x`$-axis is given by the square root of the corrected Sharvin conductance $`G_S=G_0(k_F^2\sigma ^2/4k_F𝒫/4\pi +\mathrm{\hspace{0.17em}1}/6)`$. Here $`G_0=2e^2/h`$ is the conductance quantum and $`𝒫`$ the ellipse perimeter.

The possibility to observe a stable wire in conductance histograms depends mostly on two conditions: (i) The wire has to be formed often enough to be statistically relevant; (ii) It needs to have a long enough lifetime so as to be recorded. Condition (i) depends on the ability of the system to probe new configurations, which is determined by the mobility of the atoms and the density of stable geometries in configuration space. Regarding point (ii), a calculation of the lifetime of a nanowire is beyond the scope of this paper, but one can expect that wires that are linearly stable up to larger temperatures $`T_{max}`$ will have longer lifetimes at low temperatures lif .

Based on these considerations, we extract the most stable configurations from Fig. 2, defined as the geometries that persist up to the highest temperature compared to their neighboring configurations. Those wires are marked with numbered arrows in Fig. 2. For each stability peak, we extract its mean Sharvin conductance and its width, and plot them as a function of the peak number in Fig. 3, together with the experimental data from Ref. Yanson et al. (1999). Note that the striking fit is only possible when including elliptical nanowires, for which the corresponding aspect ratios $`\epsilon `$ are also shown in Fig. 3.

The heights of the dominant stability peaks in Fig. 2 exhibit a periodic modulation, with minima occuring at peaks 4, 8, 12, 18, etc. The positions of these minima are in perfect agreement with the observed supershell structure in conductance histograms of alkali metal nanowires Yanson et al. (2000). Interestingly, the nodes of the supershell structure, where the shell effect for a cylinder is suppressed, are precisely where the most stable elliptical nanowires are predicted to occur. Thus the Jahn-Teller distortions and the supershell effect are inextricably linked.

While both the shell Yanson et al. (1999) and supershell Yanson et al. (2000) effects are accurately described by our stability analysis, our thermodynamic model does not directly address the complex dynamical process by which various contact structures form in MCBJ experiments, and therefore cannot describe every detail of the conductance histograms, such as their highly non-trivial temperature dependence.

Our stability analysis also reveals two highly deformed elliptical nanowires with conductance values of $`2G_0`$ and $`5G_0`$, marked in Fig. 2 by arrows A and B, respectively. They are expected to appear more rarely due to their reduced stability relative to the neighboring peaks, and their large aspect ratio $`\epsilon `$ (see Table 1) that renders them rather isolated in configuration space foo (b). Conductance histograms of the alkali metals often do show a conductance peak at $`5G_0`$ and a shallow shoulder at $`2G_0`$ Yanson et al. (1999, 2001). Nevertheless, these peaks were not included in the experimental analysis of shell structure in Refs. Yanson et al., 1999 and Yanson et al., 2001. Note that the definition of an experimental “peak” depends on the selection criteria and smoothing of the data.

Of the three alkali metals, Li, Na and K, Potassium has the lowest and Lithium the highest melting temperature. This suggests that at a given temperature the mobility of the atoms is highest in K and lowest in Li. This is reflected in the fact that the evidence of the highly deformed nanowires A and B is clear for K, can be seen for Na, but is barely visible for Li Yanson et al. (2001).

In addition to wire A, Fig. 2 shows a number of “superdeformed” nanowires with $`\epsilon >1.5`$. The most stable are at $`(G,\epsilon )=`$ (2, 1.65), (3, 2.25), (4, 2.8), (6, 1.65), (8, 2.0), (12, 1.65), (19, 1.55), (29, 1.55) and (41, 1.6). An aspect ratio near 2 is favorable for the shell effect Nolan and Twin (1988), as evidenced by the large gaps in the energy spectrum (Fig. 1), but the large cost in surface energy renders these wires less stable than the magic cylindrical wires and the “normal-deformed” wires listed in Table 1.

Conductance histograms have also been recorded for gold, both at low Costa-Krämer et al. (1997) and room temperatures Díaz et al. (2003); Medina et al. (2003); Mares et al. (2004), where evidence of electronic shell effects has been reported. Although the explanation of some features of Au wires, e.g. surface reconstruction, requires explicit inclusion of 5d-orbitals, the jellium model is sufficient to explain electron-shell effects. Performing a stability analysis for gold requires adjusting the deformation constraint (2) according to the surface tension for Au, ($`1.3\mathrm{N}/\mathrm{m}`$ Tyson and Miller (1977), $`\eta =0.61`$), which is higher than for Na. Therefore the results are somewhat different: Stable configurations appear at the same conductance values, but the temperature $`T_{\mathrm{max}}`$ up to which elliptical wires remain stable is larger, and these wires are less deformed than for Na (see Table 1). As a consequence, the probability to observe elliptical nanowires in experiments is enhanced. Indeed, several experimental histograms Costa-Krämer et al. (1997); Díaz et al. (2003) for Au show clear peaks at conductance values of $`2G_0`$ and $`5G_0`$.

In conclusion, we have presented a stability analysis of elliptical metal nanowires, using a jellium model, and have shown that Jahn-Teller-distorted wires can be stable. The derived sequence of stable cylindrical and elliptical geometries explains the experimentally observed shell and supershell structures for alkali metals. Deformed wires can explain additional conductance peaks observed in alkali metals and gold.

We are grateful to J. van Ruitenbeek for sharing his data on Na conductance histograms. This research has been supported by the DFG through SFB 276 and the EU Network DIENOW. CHZ and CAS acknowledge support from NSF grant DMR0312028.


text:
# Quantized charge pumping through a quantum dot by surface acoustic waves

## Abstract

We present a realization of quantized charge pumping. A lateral quantum dot is defined by metallic split gates in a GaAs/AlGaAs heterostructure. A surface acoustic wave whose wavelength is twice the dot length is used to pump single electrons through the dot at a frequency $`f=3`$GHz. The pumped current shows a regular pattern of quantization at values $`I=nef`$ over a range of gate voltage and wave amplitude settings. The observed values of $`n`$, the number of electrons transported per wave cycle, are determined by the number of electronic states in the quantum dot brought into resonance with the fermi level of the electron reservoirs during the pumping cycle.

preprint: 10/12/2003

The concept of quantized charge pumping was first suggested by Thouless Thouless83 . In an infinite one-dimensional channel with a gap in the excitation spectrum, non-interacting electrons are driven by a slowly moving periodic potential. The resulting dc current $`I`$ is quantized as $`I=ef`$, where $`e`$ is the electronic charge and $`f`$ the frequency of the moving potential. The electrons always remain in the ground-state and the mechanism is reminiscent of the Archimedean screw. In recent years, this idea has been extended to various finite systems in semiconductor nanostructures Hekking91 ; Altshuler99  and in particular charge pumps based on quantum dot (QD) devices have received considerable theoretical and experimental attention Brouwer98 ; Aleiner98 ; Shutenko00 ; Levinson01 ; Kouwenhoven91 ; Switkes99 . Accurate quantized charge pumping is of practical importance as a quantum standard for electric current Niu90 .

Quantized charge transport has been presented in metallic single electron tunneling (SET) Pothier92  and semiconductor QD Kouwenhoven91  devices where the system is in the Coulomb blockade regime. Phase-shifted AC-voltages are applied to surface gates to control the electron tunneling probability. Single electrons are transported through these devices at frequencies of roughly 10 MHz.

In this letter we report quantized charge pumping at a frequency of 3 GHz. A lateral QD is induced by metallic split gates in a GaAs/AlGaAs heterostructure and single electrons are pumped through the resonant states of this QD by the use of a surface acoustic wave (SAW). A key difference from the mechanism of the turnstile devices of Ref. Kouwenhoven91  and Ref. Pothier92  is that in our case as well as modulating the potential barriers of the QD, the traveling SAW also oscillates the electronic states within the QD relative to the Fermi level of the surrounding two-dimensional electron gas (2DEG).

This work is closely related to the single electron transport through one-dimensional channels by SAWs Shilton96 ; Talyanskii97 ; Janssen00 ; Ebbecke00 . Here, the combination of the dynamic potential due to the SAW and the static potential of the one-dimensional channel is thought to define traveling QDs, and the quantization of the current is a result of the Coulomb repulsion between electrons transported in these QDs. However, recent results have indicated that in these systems unintentional static QDs formed by impurity potentials could be involved in the transport mechanism Fletcher03 .

In Fig. 1 a scanning electron micrograph of the central part of the sample is shown together with a schematic diagram of the measurement set-up. Three independent split gates are fabricated to induce a QD of 500nm length in the 2DEG located 90 nm below the sample surface. The 2DEG is a GaAs/AlGaAs heterojunction with electron density of $`1.66\times 10^{15}`$ m-2 and mobility of 128 m2/Vs (both measured in the dark at T = 1.5 K). The two outer split gates define the entrance and the exit of the QD and are both 100 nm long with gaps of 600 nm. The middle gate is 100 nm long and has a gap of 800 nm. Two samples with identical split gate dimensions have been investigated. On sample $`A`$ there are two SAW transducers, on opposite sides of the split gates, with an interdigital period of 1 $`\mu `$m. This was chosen to give a SAW of wavelength approximately twice the QD length. The SAW velocity of 2800 m/s on GaAs gives an operating frequency of 2800 MHz. Sample $`B`$ has been processed with interdigital periods of 1 $`\mu `$m (2800 MHz) and 3 $`\mu `$m (980 MHz).

When the two outer split gates were set to a conductance $`G_{1,3}<2e^2/h`$ a series of equally spaced Coulomb blockade oscillations (CBO) was detected as the voltage on the middle gate, $`V_2`$, was swept close to conduction pinch-off. By varying the applied source-drain bias voltage, a charging energy of the induced QD of $`E_C`$0.25 meV was deduced. At the temperature of 0.3 K of our measurements no distinct structure was observed in the CBO peaks and therefore we conclude that the difference between the quantum mechanical eigenstates of the QD is $`\mathrm{\Delta }ϵ`$ 25 $`\mu `$eV $`kT`$.

Fig. 2 presents measurements of the acousto-electric current as a function of the voltage on the middle gate for sample $`A`$. The two outer split gates were held at a fixed voltage to define the entrance and exit tunnel barriers of the QD. The rf power applied to the SAW transducer was varied from $`P_{SAW}`$ = -30 dBm (top curve) to $`P_{SAW}`$ = 0 dBm (bottom curve). The potential amplitude of the SAW resulting from these drive levels has been deduced by a method described in Ref. Fletcher03 ; 0 dBm measured at the generator corresponds to an amplitude $`A_{SAW}`$ 1 mV at the QD. For the smallest SAW level ($`P_{SAW}`$ = -30 dBm, $`A_{SAW}`$ 0.03 mV), the current shows several CBO peaks before conductance-pinchoff, exactly as for no applied SAW. For larger amplitudes, the positive current due to the small applied source-drain bias changes to a negative current driven against this bias by the SAW. This acousto-electric current shows a very regular pattern which is aligned with both the CBO positions and the values of $`I=nef`$. For middle gate settings which give a peak in the CBO, the curves are most strongly grouped around odd values of $`n`$; for gate settings between peaks (maximum blockade), the curves tend to even values of $`n`$.

The piezoelectric potential of the traveling surface acoustic wave $`A(x\omega t)`$ can be regarded as a superposition of two standing waves: $`A(x\omega t)=A_0cos(2\pi t/T)sin(2\pi x/\lambda )+A_0sin(2\pi t/T+\pi )cos(2\pi x/\lambda )`$ with $`T`$ the period of the pumping cycle and $`\lambda `$ the wavelength of the SAW. The dimensions of the device have been designed such that the wavelength $`\lambda `$ is double the length of the induced QD. The resulting phase difference $`\varphi =\pi `$ in $`sin(2\pi x/\lambda )`$ and $`cos(2\pi x/\lambda )`$ at the entrance and exit potential barriers leads to a turnstile like modulation. As a result of the Coulomb blockade, a fixed number of electrons are transported during each SAW cycle. In contrast to the original turnstile mechanism presented in Ref. Kouwenhoven91  and Ref. Pothier92  our device does not require a source-drain bias to pump electrons through the QD. The SAW changes the defining potential of the static QD and therefore moves the single electron QD states energetically during the SAW cycle. This has the effect of producing a ‘local bias’ across the tunnel barriers. In particular this SAW transport mechanism allows electrons to be pumped against an applied source-drain bias and can be regarded as a nanometer scale realization of an Archimedean screw.

The periodic pattern observed in Fig. 2 can be explained by using the turnstile model and the assumption that the SAW amplitude provides a local bias. Changing the middle gate voltage brings different discrete electronic states of the QD into resonance with the Fermi levels of the source and drain 2DEG. At a CBO peak, one state is exactly aligned with the 2DEG, and the SAW transports one electron per cycle. At the gate voltages exactly in between the CBO peaks the Fermi levels of the source and drain 2DEG are in between two QD states. In this case the SAW amplitude needs to be increased to at least the level separation of the QD states. The SAW amplitude then encloses two QD states and therefore two electrons per SAW cycle are transported Fletcher03 .

In Fig. 3a) we show the acousto-electric current for a wide range of middle gate voltage, $`V_2`$, and SAW amplitude measured in sample $`A`$. Here the numerical derivative of the current with respect to applied SAW amplitude is plotted as a grey-scale map; the parameter space of Fig. 2 is outlined in the lower right corner. It can be seen that at least 18 different single-electron states of the QD have been brought into resonance by the SAW for different settings of $`V_2`$. After conduction pinch-off (the white lower left area of Fig. 3a) the pattern shifts to progressively higher SAW amplitudes. This can be explained by the fact that the conductance of the two outer point contacts decreases for more negative $`V_2`$ and hence the SAW amplitude has to be increased in order to pump electrons through the QD. The precision of the current quantization is enhanced beyond conduction pinch-off because co-tunneling and leakage as a result of the applied source-drain bias are progressively suppressed. The line graph in Fig. 3b shows a single acousto-electric current scan at a relatively high SAW amplitude to demonstrate the $`I=nef`$ sequence.

Results very similar to those shown in Fig. 3a, were obtained using the transducer on the opposite side of the QD; this transducer launches a SAW in the reverse direction, so the polarity of the pumped current is also reversed. A similar result was also obtained with sample $`B`$ using the 1 $`\mu `$m transducer. When using the second transducer on this sample to launch a SAW of 3 $`\mu `$m wavelength, the pattern of transitions remained largely unchanged; in this case, however, the acousto-electric current was less well quantized and found to be more sensitive to the applied source-drain bias. For a SAW of 3 $`\mu `$m wavelength, the phase difference between the entrance and exit barriers in a 500 nm QD is $`\pi /3`$. Therefore, during part of the pumping cycle both barriers are relatively low and the current is sensitive to co-tunneling and leakage driven by the applied source-drain bias. The periodic pattern of current quantization is similar in all cases, as this is determined by the properties of the QD, which was induced by identical split-gates on both samples.

In summary we have realized quantized charge pumping mediated by surface acoustic waves at 3 GHz which is an increase in frequency of more than 2 orders of magnitude over the previously reported turnstile devices. The use of a SAW, rather than direct gate modulation, provides a simple means to interact with a QD at such high frequencies. Although the accuracy of the current quantization observed in our first devices is only of the order $`10^3`$, this new system opens a promising route for the development of a quantum standard of current. Further experiments with smaller quantum dots should improve the accuracy. An interesting system to investigate in this context would be carbon nanotubes where a charging energy of more than 6 meV has been demonstrated Nygard01 .

###### Acknowledgements.

We thank Miriam Blaauboer, Peter Wölfle and Piet Brouwer for useful discussions. F.J.A wants to thank the Glazebrook Foundation and the Deutsche Forschungsgemeinschaft for their support. This work was supported by the National Measurement System Policy Unit of the Department of Trade and Industry, UK


text:
# The effect of trap anharmonicity on the critical temperature for Bose-Einstein condensation

## Abstract

The anharmonicity of a magnetic atomic trap at long distances from its center (forty times the radius of the first atomic orbit along the elongated axis, i.e. hundreds of microns) is shown to increase up to 45 % the temperature for Bose-Einstein condensation. This effect is perhaps small in the traditional traps, but should certainly be taken into account in magnetic microtraps, which characteristic dimensions are of order 1 mm.

Unlike the situation in liquid Helium, atomic vapors undergoing Bose-Einstein condensation (BEC) in magnetic traps are very rarefied Legget . As a consequence, the free-boson model works extremely well. Interaction effects, computed in mean-field approximation, are shown to decrease the critical temperature for BEC in a few percents Dalfovo . Still lower corrections come from the approximation of the one-particle discrete spectrum by a continuum of states, the so-called “finite-$`N`$” corrections Dalfovo , which also decrease the critical temperature.

In the present paper, we show that a relatively important increase of the critical temperature (up to 45 %) could be related to trap anharmonicities at “long” distances from the center. By long distances, we mean around forty times the radius of the first atomic orbit along the elongated direction in the trap potential, i.e. around 200 $`\mu `$m. In the commonly used traps, which characteristic dimensions are a few centimeters, the belief is that anharmonic effects should be very weak at distances of 200 $`\mu `$m. Nevertheless, we present an example of a recent experiment Pereira , in which the measured values of the total number of atoms in the trap and the critical temperature are not consistent unless anharmonicity (or other) effects raising $`T_c`$ are included. On the other hand, the trap potential should certainly be anharmonic at 200 $`\mu `$m in the recently developed atomic microtraps microtraps , which characteristic dimensions are of order 1 mm.

In the continuum free-boson model, the critical temperature is simply estimated from the conditions $`N_0=0`$, $`\mu =0`$ in:

$$NN_0=_0^{\mathrm{}}\frac{g(ϵ)\mathrm{d}ϵ}{\mathrm{exp}(\frac{ϵ\mu }{k_BT})1}.$$
(1)
$`N_0`$ is the number of atoms in the condensate, and $`ϵ0`$ is the one-atom excitation energy. The chemical potential takes values $`\mu 0`$. $`g`$ is the density of states.

Magnetic traps with cylindrical symmetry are very common. The trap potential, coming from the interaction of an hyperfine atomic species with the magnetic field, is proportional to the latter which, near the trap center, has a minimum and is written as:

$$B(\rho ,z)=B_0+b_\rho \rho ^2+b_zz^2.$$
(2)
The corresponding harmonic oscillator frequencies are $`\omega _\rho `$ and $`\omega _z`$. It is useful to define $`\omega _0=(\omega _\rho ^2\omega _z)^{1/3}`$, in terms of which the level density is written $`g(ϵ)=\frac{1}{2}ϵ^2/(\mathrm{}\omega _0)^3`$, and the critical temperature becomes:

$$k_BT_c^0=\left(\frac{N}{\zeta (3)}\right)^{1/3}\mathrm{}\omega _0=0.94N^{1/3}\mathrm{}\omega _0,$$
(3)
where $`\zeta (x)=_{n=1}^{\mathrm{}}1/n^x`$ is the Riemann Zeta function.

The corrections to $`T_c`$ due to interaction or finite-$`N`$ effects take, respectively, the form Dalfovo ; Salasnich :

$$\frac{T_cT_c^0}{T_c^0}=1.33\frac{a}{r_0}N^{1/6},$$
(4)
$$\frac{T_cT_c^0}{T_c^0}=0.73\frac{\overline{\omega }}{\omega _0}N^{1/3},$$
(5)
where $`r_0=\sqrt{\mathrm{}/(m\omega _0)}`$, $`a`$ is the s-wave scattering length, and $`\overline{\omega }=(\omega _z+2\omega _\rho )/3`$. Typical values $`N10^6`$, $`a/r_010^2`$ lead to small numbers in the r.h.s. of Eqs. (4,5). As mentioned, the corrections (4) are computed within mean-field theory. More elaborated variational calculations based on diffusion Monte-Carlo Giorgini  or hypernetted chain theory Fantoni  support the validity of mean-field approximations at the relevant particle densities. Notice that both corrections (4) and (5) are small and make $`T_c`$ decrease. In a general power-law trap, the corrections may be of any sign and their relative magnitudes may be still greater Salasnich . However, for usual magnetic traps one expects a harmonic potential coming from (2), at least in the vicinity of the trap center.

Let us suppose that at long distances, $`rr_0`$, anharmonic terms are relevant in (2). In fact, one expects a saturation of the magnitude of $`B`$, which should approach the bias field value as the distance increases. In a trap potential that saturates at long distances, the density of energy levels decreases at high energiesnota . It means that the occupation of the ground state is increased, and thus the critical temperature is raised.

We will model the density of levels of the (saturating) anharmonic potential in the simplest way: by truncating the one-particle spectrum. It means that an upper integration limit, $`ϵ_{max}`$, is set in (1). The critical temperature is thus determined from:

$`t^3\zeta (3)`$ $`=`$ $`\zeta (3)Li_3(e^{tx})txLi_2(e^{tx})`$ (6)
$`+`$ $`{\displaystyle \frac{(tx)^2}{2}}\mathrm{ln}(1e^{tx}),`$
where we have defined $`t=T_c^0/T_c`$, $`x=ϵ_{max}/(k_BT_c^0)`$, and $`Li_k(z)=_{n=1}^{\mathrm{}}z^n/n^k`$ is the polylogarithm function.

The corresponding correction to $`T_c^0`$, i. e. $`1/t1`$, is shown in Fig. 1 as a function of $`x`$. It is around 10% when $`ϵ_{max}4k_BT_c^0`$. As $`k_BT_c^010^3\mathrm{}\omega _z`$ in many of the experiments, energies of order $`4k_BT_c^0`$ involve particle orbits of radius around 60 $`r_{0z}`$. On the other hand, when $`ϵ_{max}/k_B2.5T_c^0`$, corrections to $`T_c^0`$ reach 40 %, and the cutoff distance is around 40 $`r_{0z}`$. A fit to the long-x tail of the curve leads to:

$$\frac{T_cT_c^0}{T_c^0}=22\left(\frac{k_BT_c^0}{ϵ_{max}}\right)^4.$$
(7)
In quality of example, let us consider the trap used in Ref. Pereira, , in which BEC of spin-polarized Helium was achieved. The example is quite interesting because independent measurements of $`N`$ and $`T_c`$ are reported. The authors-provided frequencies are $`\nu _z=115`$ Hz, $`\nu _\rho =1090`$ Hz. The critical temperature, number of atoms in the trap and the scattering length are estimated as $`T_c=4.7\pm 0.5\mu `$K, $`N=(5\pm 2.5)\times 10^6`$, $`a16`$ nm, respectively. Turning back to Eqs. (3-5), we obtain that a temperature $`T_c=4.7\mu `$K, corresponds to a number of atoms $`N=1.4\times 10^7`$, well above the error bars. On the other hand, if we take the lowest value $`T_c=4.2\mu `$K, then the number of atoms is $`N=8\times 10^6`$, a value closer but still outside error bars. The reason for this apparent discrepancy may be either the rough estimation of $`N`$, or the anharmonicity corrections, which have not been included.

In conclusion, we have shown that the saturation of the magnetic field at distances around 200 - 300 $`\mu `$m from the trap center leads to a decrease of the density of energy levels, and thus to an increase of the ground-state occupation and to an increase of $`T_c`$. We have modeled the decrease of the level density simply by means of an energy cutoff in the single-particle spectrum. The corresponding correction to $`T_c^0`$ partially cancels out with interaction or finite-$`N`$ effects, both of which lower the critical temperature. Although expected to be small, detailed calculations of the magnetic field configurations in the commonly used traps are needed in order to evaluate the role of anharmonicity corrections. With regard to the atomic microtraps, which characteristic dimensions are of order 1 mm, anharmonicities are expected to be very important and should be taken into account.

###### Acknowledgements.

Part of this work was carried out at the Abdus Salam ICTP. The author acknowledges the ICTP Associate and Federation Schemes for support.


text:
# A theory of growth by differential sedimentation, with application to snowflake formation

## I Introduction

Simple models of cluster-cluster aggregation have been the focus of a great deal of interest, particularly over the last two decades. The structure of aggregates formed through a variety of dominating mechanisms (eg. diffusion limited meakin , reaction limited balletal  and ballistic motion jullienandkolb ) have been studied through theoretical, experimental, and computational work.

Another aggregation mechanism which is relevant to several physical systems is that of differential sedimentation. Particles with a range of size and/or shape will almost inevitably sediment through a fluid at different speeds under the influence of gravity, leading to collisions. If there is some mechanism by which the particles stick on contact then aggregates will be formed. An example of this kind of phenomenon is the aggregation of ice crystals in Cirrus clouds. Small ‘pristine’ ice particles are formed at the top of the cloud, and proceed to fall through it, colliding with one another and sticking to produce aggregates (snowflakes).

The aim of this paper is to provide a simple model for growth by differential sedimentation which captures the essential physics of the system in the inertial flow regime, and to consider its application to snowflake formation. It is divided into five main parts - a description of the model and the assumptions underlying it; details of computer simulations and the results obtained from them; a theory section which offers an argument to account for the behaviour observed in the simulations; an investigation of the model’s applicability to snowflake formation; and a concluding discussion. The simulation results and their comparison to real cloud data have been presented briefly in a separate letter ourshortpreprint .

## II Model

We focus on the dilute limit, where the mean free path between cluster-cluster collisions is large compared to the nearest neighbour distance between clusters. In this regime we can limit our interest to individual binary collision events, ignoring spatial correlation. As further simplifying approximations, we assume that clusters have random orientations which do not significantly change during a close encounter, that collsion trajectories are undeflected by hydrodynamic interaction, and that any cluster-cluster contacts result in a permanent and rigid junction.

In order to sample the collisions between clusters, we first formulate a rate of close approach. For any two clusters $`i`$, $`j`$ with nominal radii (see below) $`r_i`$ and $`r_j`$ respectively and fall speeds $`v_i`$ , $`v_j`$ , the frequency with which their centres pass closer than a distance $`(r_i+r_j)`$ is proportional to the total area over which trajectories yielding a close approach event are possible, and the relative speed of the pair. This is illustrated in figure 1. The rate constant for approach closer than centre-to-centre separation $`(r_i+r_j)`$ is therefore given by:

$$\mathrm{\Gamma }_{ij}=\pi (r_i+r_j)^2\left|v_iv_j\right|.$$
(1)
In our computer simulations the nominal radii are chosen to fully enclose each cluster and the close approach rate calculated above is exploited to preselect candidate collision events. Collisions are accurately sampled by sequentially choosing pairs of clusters with probability proportional to $`\mathrm{\Gamma }_{ij}`$, checking each pair for collision along one randomly sampled close approach trajectory, and correspondingly joining that cluster pair if they do indeed collide. In the theoretical arguments presented in section four, we make the simplifying assumption that all close approaches lead to collisions (or at least a fixed fraction of them do), using nominal radii based on fractal scaling from the cluster masses.

The model is completed by an explicit form for the fall speeds entering equation (1). We assume that the clusters are at most only partially penetrated by the fluid flow past them, so that cluster radius is the relevant length governing the drag force law. Then qualitatively and by dimensional argument we expect the same drag behaviour as for a falling sphere, which may be written in the form:

$$F_d=\rho \nu _kf(R_e)$$
(2)
where $`f`$ is a function of the Reynolds number $`R_e=rv/\nu _k`$ alone, $`\rho `$ is the density of the surrounding fluid, and $`\nu _k`$ is the kinematic viscosity. Although details of the function $`f(R_e)`$ should be different from spheres, we still expect to have inertial and Stokes regimes where $`f`$ takes the forms:

$`f(R_e)\left\{\begin{array}{c}R_e^2\text{ for inertial flow}\\ R_e\text{ for viscous flow}\end{array}\right\}.`$ (5)
We consider below the general form $`f(R_e)R_e^{1/\alpha }`$, with $`\alpha `$ as an adjustable parameter in order gain understanding spanning the two extreme regimes. Setting the drag force equal to the weight $`mg`$ of the cluster, the terminal velocity is then given by

$$v\frac{\nu _k}{r}\left(\frac{mg}{\rho \nu _k^2}\right)^\alpha $$
(6)
where $`\alpha =\frac{1}{2}`$ for inertial flow and $`\alpha \text{=1}`$ for viscous flow. A more complete discussion of the fall speed is given by Mitchell mitchell96 , but provided $`d_f2`$ so that cluster projected area scales as the square of cluster radius, this reduces to a simple crossover between the above limits. The empirical crossover is very slow, spread over some three decades of Reynolds number, so fixed intermediate values of $`\alpha `$ can reasonably approximate behaviour over a significant range mitchell96 . In our simulations we took the radius determining the fall velocity to be proportional to the radius of gyration, and in our theoretical calculations we simply used the same nominal radii as for the collision cross sections above.

## III Computer simulations and results

The primary particles used at the beginning of the simulations were rods of zero thickness, half of which had a length (and mass) of unity, and half of which were twice as long and massive. Purely monodisperse initial conditions are not possible in this model, since $`|v_iv_j|`$ would be zero. Apart from this special case however, it is anticipated that the asymptotic behaviour of the system should be insensitive to the initial distribution, and indeed the results described in this section appear to be preserved for a variety of starting conditions.

In aggregation models it is typically the case (eg. Vicsek and Family vicsekandfamily ) that after the distribution has had time to ‘forget’ its initial conditions it will approach a universal shape. This is usually expressed by the ‘dynamical scaling’ ansatz, which states that as $`m,s\mathrm{}`$:

$$n_m(t)=s(t)^2\varphi \left[\frac{m}{s(t)}\right]$$
(7)
where $`n_m(t)`$ is the number of clusters of mass $`m`$ at time $`t`$, and the rescaled distribution $`\varphi `$ is a function of $`x=m/s(t)`$ alone. The quantity $`s(t)`$ is a characteristic cluster mass, and for non-gelling systems one expects that a suitable choice is given by the weight average cluster mass, $`s(t)=_im_i^2/_im_i`$. Using this choice our simulation data conform well to scaling, as shown in the left panel of figure 2.

The shape of the rescaled distribution was studied. A plot of $`_x^{\mathrm{}}\varphi (x^{})\text{d}x^{}`$ as a function of $`x`$ is shown in the right panel of figure 2 and shows an exponential decay for very large $`x`$, with a ‘super-exponential’ behaviour taking over as $`x`$ approaches unity from above. This behaviour appears to be universal for all values of $`\alpha `$ in the range studied.

For $`x1`$ the qualitative form of $`\varphi (x)`$ was found to fall into two distinct catagories depending on the value of $`\alpha `$. For $`\alpha \frac{1}{2}`$ the distribution appears to diverge as a power law: $`\varphi (x0)x^\tau `$, as shown in figure 2 for $`\alpha =0.55`$. The exponent $`\tau `$ was found to be approximately constant at $`\tau 1.6\pm 0.1`$ over the range $`\frac{1}{2}\alpha \frac{2}{3}`$. For $`\alpha <\frac{1}{2}`$ the distribution was found to be peaked, with a maximum at some small size $`x_m`$, followed by a power law decay for $`x_mx1`$.

Comparison with other aggregation models suggests that the clusters produced are likely to be fractal in their geometry, and in particular cluster mass and (average) radius should be in a power law relationship $`mr^{d_f}`$ where $`d_f`$ is the fractal dimension. A log plot of radius of gyration against mass for all the clusters produced over the course of the simulation is shown in figure 3. Also shown in this figure is the logarithmic derivative of the above plot, which shows the variation in the apparent fractal dimension of the clusters with size. From this plot, it seems that the fractal dimension approaches an asymptotic value as $`m\mathrm{}`$; in the case shown ($`\alpha =0.55`$) we estimate this value as $`d_f2.2\pm 0.1`$. The value of this limiting fractal dimension was found to vary with $`\alpha `$ as shown in figure 4. Note that our assumption $`d_f2`$, required to support the assumed fall speed relationship, is indeed satisfied for the physical range $`\alpha 1/2`$.

## IV Theory

The most common theory used to describe cluster-cluster aggregation problems is that of von Smoluchowksi smoluchowski , which provides a set of mean-field rate equations for the evolution of the cluster mass distribution:

$$\frac{\text{d}n_k(t)}{\text{d}t}=\frac{1}{2}\underset{i+j=k}{}K_{ij}n_i(t)n_j(t)n_k(t)\underset{j=1}{\overset{\mathrm{}}{}}K_{ij}n_j(t).$$
(8)
where $`n_k(t)`$ is the number of clusters of mass $`k`$ at time $`t`$ (per unit volume). The kernel $`K_{ij}`$ contains the physics of the problem, being a symmetric matrix, the elements of which govern the rate of aggregation between pairs of clusters expressed (only) in terms of their masses $`i`$ and $`j`$. Analytical solutions of Smoluchowski’s equations have not been obtained except for a few special cases of $`K_{ij}`$. However, Van Dongen and Ernst vdande85  have shown that for non-gelling kernels (see below) the solutions approach the dynamical scaling form of equation (7) in the large-mass, large-time limit; substituting this into equations (8) allows one to obtain some information about the asymptotic behaviour of the rescaled cluster size distribution $`\varphi (x)`$.

To apply this theory we need to compute the reaction rates $`K_{ij}`$ , which means averaging collision rate with respect to cluster geometry at fixed masses. This we estimate by substituting averages from fractal scaling for the radii in equations (1) for the close approach rate and (6) for the fall speeds, and assuming constant collision efficiency leading to:

$$K_{ij}\left|i^{\alpha 1/d_f}j^{\alpha 1/d_f}\right|\left(i^{1/d_f}+j^{1/d_f}\right)^2.$$
(9)
Van Dongen and Ernst’s analysis is sensitive to two exponents characterising the scaling of the coagulation kernel in the limit $`1ij`$,

$$K_{ij}i^\mu j^\nu $$
(10)
which in our case yields:

$`\mu `$ $`=`$ $`\mathrm{min}(0,\alpha d_f^1)`$ (11)
$`\nu `$ $`=`$ $`\mathrm{max}(\alpha +d_f^1,2d_f^1).`$ (12)
A third exponent combination $`\lambda =\mu +\nu =\alpha +d_f^1`$ controls the growth of the average cluster mass through the differential equation $`\dot{s}(t)=ws(t)^\lambda `$, where $`w`$ is a constant, and for the non-gelling case we require $`\lambda 1`$.

Our identification of the exponent $`\nu `$ is crucial to a mechanism by which the fractal dimension can control the dynamics. If the fractal dimension is low enough, then the exponent $`\nu `$ will exceed unity. However, Van Dongen vandongen  has shown that the Smoluchowski equations predict the formation of an infinite cluster instantly in such a situation. In a finite system this clearly cannot occur, and it simply means that a few clusters will quickly become much larger than the others with their growth dominated by accretion of small ones. In this scenario the growth of the large clusters approaches that of ballistic particle-cluster aggregation, where it has been shown by Ball and Witten ballandwitten  that the fractal dimension of the clusters produced tends to $`d_f=3`$. This increased fractal dimension reduces the value of $`\nu `$, forcing it back to a value of one if $`\alpha \frac{2}{3}`$. Through this feedback mechanism, a bound is placed on the fractal dimension $`d_f\mathrm{max}[2,(1\alpha )^1]`$ for $`\alpha \frac{2}{3}`$.

The system could perhaps settle in a state where $`\nu <1`$. However, the growth in such a regime is much less biased towards collisions between clusters of disparate sizes, and the distribution is relatively monodisperse. This would tend to make collisions between clusters of a similar size likely, leading to much more open structures, with a lower fractal dimension, in turn acting as a feedback mechanism to increase the value of $`\nu `$. The authors suggest that, at least over some range of $`\alpha `$, this effect will force the system towards the $`\nu =1`$ state. The discontinuity in the polydispersity of the system at $`\nu =1`$ forces the system to organise itself such that it can remain at that point. This is similar to the argument put forward by Ball et al balletal  for reaction limited aggregation.

If it is accepted that $`\nu 1`$ then the fractal dimension of the clusters produced ought to be directly predictable from equation (12) :

$$d_f=\mathrm{max}[2,(1\alpha )^1],\alpha \frac{2}{3}.$$
(13)
A curve showing this theoretical behaviour is superimposed on the simulation data in figure 4, and appears to show good agreement up to $`\alpha \frac{2}{3}`$. For $`\alpha >\frac{2}{3}`$ the theoretical prediction is that $`d_f=3`$ and $`\nu =\alpha +\frac{1}{3}>1`$, but because of its somewhat pathological nature we have not attempted to make simulations in this regime. It is however clear from the extrapolation of our results in figure 4 that this is likely to hold.

Obtaining an exact form for the cluster size distribution $`\varphi (x)`$ is a non-trivial exercise. However, following the methodology of Van Dongen and Ernst vdande85 , we consider the small-$`x`$ behaviour of $`\varphi (x)`$ when $`d_f<\alpha ^1`$ (ie. $`\mu <0`$). In such a regime the small-$`x`$ behaviour is dominated by collisions between clusters of disparate sizes; the gain term in the Smoluchowski equations may therefore be neglected, and one attempts to solve the integro-differential equation: $`w[x\varphi ^{}(x)+2\varphi (x)]=\varphi (x)_0^{\mathrm{}}K(x,y)\varphi (y)\text{d}y`$. For $`xy`$, the kernel (9) may be approximated to $`K(x,y)x^\mu y^\nu y^\lambda `$, and one obtains:

$$\varphi (x)=x^\tau \mathrm{exp}\left[\frac{x^\mu p_\nu }{w\mu }\right]$$
(14)
where $`p_i`$ is the $`i^{th}`$ moment of the rescaled distribution $`\varphi (x)`$, and the exponent $`\tau `$ is given by $`\tau =2+p_\lambda w^1`$. It is clear that $`lim_{x0}[\varphi (x)]=0`$. As $`x`$ increases from zero, $`\varphi (x)`$ also increases, until reaching a maximum at $`x_m=(w\tau /p_\nu )^{1/\mu }`$. For $`x_mx1`$ the distribution has an approximately algebraic decay $`\varphi (x)x^\tau `$. This ‘bell-shaped’ curve is consistent with the behaviour seen in the computer simulations when $`\alpha <\frac{1}{2}`$.

In the case $`d_f>\alpha ^1`$, it has been shown vdande85  that for all kernels with $`\mu =0`$, $`\nu 1`$ the cluster size distribution diverges as $`x0`$ with the form

$$\varphi (x)x^\tau $$
(15)
where $`\tau =2p_\lambda w^1`$. This behaviour is consistent with the simulation for $`\alpha \frac{1}{2}`$. The change in the qualitative shape of $`\varphi (x1)`$ around $`\alpha =\frac{1}{2}`$ then is further evidence to suggest that the system selects to sit at $`\nu =1`$.

The shape of $`\varphi (x1)`$ has also been studied by Van Dongen and Ernst vdande87 . They have shown that for non-gelling kernels, the tail of the distribution is expected to take the form

$$\varphi (x)x^\theta \text{e}^{\delta x}$$
(16)
where $`\theta `$ and $`\delta `$ are constants. This would appear to be consistent with the behaviour observed in the simulations for all values of $`\alpha `$, providing an exponentially dominated cut-off at large $`x`$.

## V Application to snowflake formation

The principle motivation for the model presented in this paper was to attempt to understand some of the properties of Cirrus clouds. These are high altitude clouds with a base betwen 5,500 and 14,000 metres and they are usually composed solely of ice crystals PandKlett . Amongst others, Heymsfield and Platt heymandplatt  have observed that the ice crystals in these clouds are predominantly composed of columns, bullets, bullet-rosettes and aggregates of these crystal types. It is these aggregates which we hope to model, since the dominant mechanism by which they grow is believed to be through differential sedimentation (eg. Field and Heymsfield field ). We therefore ignore the effects of diffusional growth, turbulence, mixing, and particle breakup, in order to concentrate on the effects of this mechanism alone. The Reynolds number for aggregates of a few crystals is typically between $`10`$$`100`$ which ought to be modelled acceptably by our inertial flow approximation. Because we have not modelled the detailed hydrodynamics we may also be ignoring subtleties such as wake capture.

All of the results below are presented for the purely inertial regime (assumed to be the most relevant to this problem) where $`\alpha =\frac{1}{2}`$. The initial particles were rods of zero thickness - however, the asymptotic behaviour is anticipated to be insensitive to the initial conditions, and indeed by running the simulation with ‘bullet rosettes’ for the initial particles (three rods, crossing one another at right angles, through a common centre), no change in the end results were found, only in the approach to scaling.

Ice crystal aggregates have been studied through the use of cloud particle imagers during aircraft flights through ice clouds. Sample images from such a flight are shown in figure 5, alongside some of our simulation clusters. Using this experimental data, the geometry and size distribution of ice particles in these clouds has been studied, allowing for quantitative comparison between theory and experiment.

The fractal dimension of snowflakes in Cirrus may be inferred from the work of Heymsfield et al heyma . By measuring the effective density $`\rho _e`$ of bullet and bullet-rosette aggregates as a function of their maximum linear dimension $`D`$, and fitting a power law to their data, they found the relationship $`\rho _eD^{0.96}`$. This scaling implies that the aggregates have a fractal dimension of $`d_f=2.04`$, which is consistent with the values predicted by our model (simulation giving $`d_f=2.05\pm 0.1`$ and theory giving $`d_f=2`$).

The aspect ratio of the clusters may also be calculated. Random projections of simulation clusters were taken. The maximum dimension of the projection $`D`$ was measured, as was the maximum dimension in the direction perpendicular to that longest axis, $`D_w`$. The ratio of these two spans were binned by maximum dimension, averaged, and plotted as a function of $`D`$ as shown in figure 6. The ratio quickly approaches an asymptotic value of approximately $`0.65\pm 0.05`$. This compares well to the measurements of Korolev and Isaac korolevandisaac , where the ratio seems to approach a value of $`0.60.7`$.

Finally, the shape of the snowflake distribution of linear size may also be compared with experiment. Field and Heymsfield field  presented particle size distributions of the maximum length $`D`$ of ice particles in a Cirrus cloud. The data were obtained with an aircraft and represent in-cloud averages of particle size distributions (number per unit volume per particle size bin width) along 15km flight tracks ranging from an altitude of 9500m ($`50^{}`$C) to 6600m ($`28^{}`$C). To compare this data to the distributions obtained from simulation, we first normalise the data, and then make use of the dynamical scaling form (7), to collapse the distributions onto a single curve. Details of this are given in the appendix to this paper. The resulting histograms are shown in figure 7 and appear to show quite good agreement, given the level of approximation present in our model.

## VI Discussion and conclusions

A simple mean-field model of aggregation by differential sedimentation of particles in an inertial flow regime has been constructed, simulated by computer, and analysed theoretically in terms of the Smoluchowski equations. It has been shown that there is strong numerical evidence, in addition to a theoretical argument, to back up the idea that the polydispersity of the distribution and the fractal dimension feed back on one another in such a way as to stabilise the system at $`\nu =1`$. Above this value, the dominance of collisions between clusters of very different sizes is so great as to push $`d_f`$ towards a value of three. This in turn pulls the exponent $`\nu `$ back down to unity. For $`\nu <1`$ the system is quite monodisperse, resulting in relatively many collisions between clusters of similar sizes, and the fractal dimension is reduced, forcing $`\nu `$ back up. The discontinuity in the shape of the distribution around $`\nu =1`$ is thought to provide the mechanism by which the system can stabilise at that point.

If it is accepted that $`\nu 1`$, then the fractal dimension of the clusters produced may be predicted, and figure 4 shows that this prediction agrees well with simulation results for $`0\alpha \frac{2}{3}`$. The sudden change in the behaviour of $`d_f(\alpha )`$ and in the small-$`x`$ form of the cluster size distribution around $`\alpha =\frac{1}{2}`$ is strong evidence for the self-organisation proposed between $`d_f`$ and $`\nu `$.

For $`\alpha >\frac{2}{3}`$ the system is forced into a regime where $`\nu >1`$, which has been regarded as unphysical because the Smoluchowski equation (8) predicts infinite clusters in zero time vandongen . In the light of our results this regime merits further study beyond the Smoluchowski equation approximation workinprogress . The value $`\alpha =1`$ is given by viscous flow, but here our form for $`\mathrm{\Gamma }_{ij}`$ does not include all of the relevant physics: in particular, small clusters may be caught in the fluid flow, and swept around larger clusters rather than hitting them, reducing the dominance of big-little collisions. This has been discussed in more detail for the particle-cluster aggregation case by Warren et al warrenetal .

The application of the model to the formation of ice crystal aggregates in Cirrus clouds has been considered: the fractal dimension, aspect ratio, and shape of the cluster size distribution seen in the model were all found to be consistent with experimental data. This is a promising indication that the ideas presented in this paper may be an acceptable model for the essential physics of snowflake aggregation in Cirrus.

###### Acknowledgements.

We are grateful to Aaron Bansemer for initial processing of the size distribution data and Carl Schmitt for supplying the CPI images. This work was supported financially by the Engineering and Physical Sciences Research Council, and The Meteorological Office.

*

## Appendix A Scaling of the cluster radius distribution

Experiments have reported the distribution of ice aggregates by linear span rather than by mass, and we present here how that distribution $`\frac{\text{d}N}{\text{d}D}`$ should naturally be rescaled. This tests the dynamical scaling ansatz which, for the mass distribution, gave $`\frac{\text{d}N}{\text{d}m}=n_m=s^a\varphi (m/s)`$, where $`a=2`$ in mass-conserving systems. We anticipate fractal scaling so that $`mD^{d_f}`$ and hence:

$$\frac{\text{d}N}{\text{d}D}D^{d_f1}s^a\varphi \left(\frac{m}{s}\right).$$
(17)
From this expression we may calculate the moments of the distribution $`M(b)\frac{\text{d}N}{\text{d}D}D^b\text{d}D`$ in terms of the average cluster mass $`s(t)`$:

$$M(b)s^{a+1+b/d_f}_{1/s}^{\mathrm{}}x^{b/d_f}\varphi (x)\text{d}x$$
(18)
where $`x=m/s`$. At small sizes we expect $`\varphi (x)x^\tau `$. If $`b>d_f(\tau 1)`$ therefore, the integral converges as $`s\mathrm{}`$, and $`M(b)s^{a+1+b/d_f}`$. From our simulations, we have measured $`\tau 1.6`$, $`d_f2`$, and so the lowest integer moment which scales in this way is the second. We therefore choose this to normalise our data:

$$[M(2)]^1\frac{\text{d}N}{\text{d}D}D^{d_f1}s^{12/d_f}\varphi \left(\frac{m}{s}\right)$$
(19)
which, defining the average cluster diameter $`D^{}M(3)/M(2)s^{1/d_f}`$ yields:

$$[M(2)]^1\frac{\text{d}N}{\text{d}D}(D^{})^3\psi \left(\frac{D}{D^{}}\right),$$
(20)
where $`\psi (y)=y^{d_f1}\varphi (y^{d_f})`$. Hence, if we assume that $`d_f`$ approaches a constant value, plots of $`\{[M(2)]^1.\frac{\text{d}N}{\text{d}D}.(D^{})^3\}`$ against $`(D/D^{})`$ should all lie on a single curve.


text:
# Sagnac interferometry based on ultra-slow polaritons in cold atomic vapors

## Abstract

The advantages of light and matter-wave Sagnac interferometers – large area on one hand and high rotational sensitivity per unit area on the other – can be combined utilizing ultra-slow light in cold atomic gases. While a group-velocity reduction alone does not affect the Sagnac phase shift, the associated momentum transfer from light to atoms generates a coherent matter-wave component which gives rise to a substantially enhanced rotational signal. It is shown that matter-wave sensitivity in a large-area interferometer can be achieved if an optically dense vapor at sub-recoil temperatures is used. Already a noticeable enhancement of the Sagnac phase shift is possible however with much less cooling requirements.

The relative phase accumulated by two counter-propagating waves in a rotating ring interferometer, first observed by Sagnac Sagnac-1913 , is a well established tool to detect intrinsically the rotation of a system. Two types of Sagnac gyroscopes have been developed which have very different characteristics depending on the type of wave phenomena they are based on Stedman-RPP-1997 . Optical devices achieve high rotational sensitivity because of the large interferometer area Chow-RMP-1985 . On the other hand, the rotational sensitivity per unit area of matter-wave gyroscopes Colella-PRL-1975 ; Riehle-PRL-1991  exceed that of optical ones by the ratio $`mc^2/\mathrm{}\omega 10^{11}`$. They suffer, however, from the smallness of the achievable loop areas and only recently short-time rotational sensitivities have been demonstrated for matter-wave interferometers which are comparable to state-of-the-art laser gyroscopes Gustavson-PRL-1997 ; McGuirk-PRL-2000 . We here show that it should be possible to combine the large rotational sensitivity per unit area of matter-waves and the large area typical for optical gyroscopes utilizing the coherence and momentum transfer associated with ultra-slow light in cold atomic vapors with electromagentically induced transparency (EIT) Harris-PhToday-1997 ; slow-light . As the slow-down of light in EIT media is based on the rotation of dark-state polaritons from electromagnetic to atomic excitations Fleischhauer-PRL-2000 , light waves can coherently be transformed into matter-waves, if the excitation transfer is accompanied by momentum transfer and the atoms are allowed to move freely. Since the transfer is coherent and reversible it serves as a basis for a hybrid light-matter-wave Sagnac interferometer.

The rotational phase shift in an optical Sagnac interferometer is given by

$`\mathrm{\Delta }\varphi _{\mathrm{light}}={\displaystyle \frac{4\pi }{\lambda c}}𝛀𝐀,`$ (1)
where $`𝛀`$ and $`𝐀`$ are the vectors of the angular velocity and the loop area. Discussing Fresnel dragging in EIT media Leonhardt and Piwnicki suggested in Leonhardt-PRA-2000  that the reduction of the group velocity in a solid attached to the rotating body should lead to the replacement of the vacuum speed of light $`c`$ in (1) by the group velocity $`v_{\mathrm{gr}}`$ and thus to an enormous increase of rotational sensitivity. However, as will be shown here, the Sagnac phase shift in a solid medium is always given by the vacuum expression irrespective of the group velocity. In fact in the history of Sagnac interferometry there had been a longer discussion about the effect of refractive materials in the beam path of passive optical devices. It was finally recognized that materials that change the phase velocity of light have no effect on the Sagnac phase shift (see e.g. Stedman-RPP-1997 ). However, if the coherence transfer from light to medium, associated with any group-velocity reduction, is accompanied by a momentum transfer to freely moving particles a traveling matter-wave component is created which can lead to a substantial enhancement of the rotational sensitivity.

Let us consider a circular light interferometer of radius $`R`$ with an atomic vapor cell or trap in the beam path as shown in Fig. 1 attached to a rotating body. The product of angular velocity $`\mathrm{\Omega }`$ and radius $`R`$ is assumed to be sufficiently small compared to $`c`$ such that non-relativistic quantum mechanics applies Anandan-PRD-1981 . Light propagation as well as center-of-mass motion of the atoms shall be confined to the periphery of the loop, e.g. by means of optical fibers and appropriate atom potentials. Distances along the periphery will be denoted by the coordinate $`z`$. As indicated in Fig. 1 the atoms shall have three internal states $`|1,|2,|3`$ allowing for a $`\mathrm{\Lambda }`$-type Raman coupling of the $`|1|2`$ transition to the probe field described by the Rabi-frequency $`\mathrm{\Omega }_p(z,t)`$ and of the $`|2|3`$ transition to the control field, characterized by the Rabi-frequency $`\mathrm{\Omega }_c`$. For the present application $`\mathrm{\Omega }_c`$ is assumed to be much larger than $`\mathrm{\Omega }_p`$.

Under two-photon resonance conditions, the control field generates electromagnetically induced transparency for the probe field associated with a reduction of the group velocity according to $`v_{\mathrm{gr}}=c\mathrm{cos}^2\theta `$ with $`\mathrm{tan}^2\theta g^2n/|\mathrm{\Omega }_c|^2`$ where $`n`$ is the density of atoms and $`g=d\sqrt{\omega _p/(2\mathrm{}ϵ_0F)}`$ is a coupling constant containing the dipole moment $`d`$ of the probe transition and the transversal cross section $`F`$ of the probe beam. As has been demonstrated in several experiments slow-light  rather substantial group-velocity reductions can be achieved.

The atoms are described by the matter fields $`\mathrm{\Psi }_1(z,t)`$, $`\mathrm{\Psi }_2(z,t)`$, and $`\mathrm{\Psi }_3(z,t)`$ corresponding to the three internal states. Since the sources of the fields as well as the experimental apparatus are attached to the rotating frame, the dynamics of the probe light $`\mathrm{\Omega }_p`$ and the matter-wave components $`\mathrm{\Psi }_\mu `$ will be described in the rotating frame. For this the total Hamiltonian in a non-rotating frame $`H`$ will be transformed according to $`\mathrm{exp}\{iGt\}H\mathrm{exp}\{iGt\}`$, where $`G=\mathrm{\Omega }\widehat{M}`$ is the generator of a uniform rotation Bialynicki-Birula-PRL-1997 . $`\widehat{M}`$ is the total (orbital) angular momentum of matter fields and light in the direction of $`𝛀`$. This leads to the following equations of motion

$`𝒟\mathrm{\Psi }_1`$ $`=`$ $`\mathrm{}\mathrm{\Omega }_p^{}\mathrm{e}^{i(\omega _ptk_pz)}\mathrm{\Psi }_2,`$ (2)
$`𝒟\mathrm{\Psi }_2`$ $`=`$ $`\mathrm{}(\omega _2i\gamma )\mathrm{\Psi }_2+\mathrm{}\mathrm{\Omega }_p\mathrm{e}^{i(\omega _ptk_pz)}\mathrm{\Psi }_1`$
$`+\mathrm{}\mathrm{\Omega }_c\mathrm{e}^{i(\omega _ctk_c^{}z)}\mathrm{\Psi }_3,`$
$`𝒟\mathrm{\Psi }_3`$ $`=`$ $`\mathrm{}\omega _3\mathrm{\Psi }_3+\mathrm{}\mathrm{\Omega }_c^{}\mathrm{e}^{i(\omega _ctk_c^{}z)}\mathrm{\Psi }_2,`$ (4)
with

$`𝒟`$ $``$ $`i\mathrm{}_t+{\displaystyle \frac{\mathrm{}^2}{2m}}_z^2i\mathrm{}\mathrm{\Omega }R_z.`$
Here $`\omega _p`$ and $`\omega _c`$ are the probe and coupling frequencies in the rotating frame, which is the laboratory frame, and $`k_p=\omega _p/c`$ and $`k_c^{}`$ the corresponding wave-vector components in $`z`$ direction. $`\gamma `$ describes losses by spontaneous emission from the excited state $`|2`$. We assume that initially, i.e. without applying the probe field all atoms are in state $`|1`$.

We proceed by introducing slowly-varying matter-wave amplitudes $`\mathrm{\Psi }_1=\mathrm{\Phi }_1`$, $`\mathrm{\Psi }_2=\mathrm{\Phi }_2\mathrm{e}^{i(\omega _ptk_pz)}`$, and $`\mathrm{\Psi }_3=\mathrm{\Phi }_3\mathrm{e}^{i((\omega _p\omega _c)t\mathrm{\Delta }kz)}`$, with $`\mathrm{\Delta }k=k_pk_c^{}`$ and assume single- and two-photon resonance of the carrier frequencies, i.e. $`\omega _p=\omega _2+\mathrm{}k_p^2/2m`$ and $`\omega _p\omega _c=\omega _3+\mathrm{}(\mathrm{\Delta }k)^2/2m`$. Neglecting second-order derivatives within a slowly-varying envelope approximation we arrive at the set of equations

$`\left(_t\mathrm{\Omega }R_z\right)\mathrm{\Phi }_1`$ $`=`$ $`i\mathrm{\Omega }_p^{}\mathrm{\Phi }_2,`$ (5)
$`\left(_t(\mathrm{\Omega }Rv_{\mathrm{rec}})_z\right)\mathrm{\Phi }_2`$ $`=`$ $`+(ik_p\mathrm{\Omega }R\gamma )\mathrm{\Phi }_2`$ (6)
$`i\mathrm{\Omega }_c\mathrm{\Phi }_3i\mathrm{\Omega }_p\mathrm{\Phi }_1,`$
$`\left(_t(\mathrm{\Omega }R\eta v_{\mathrm{rec}})_z\right)\mathrm{\Phi }_3`$ $`=`$ $`+i\eta k_p\mathrm{\Omega }R\mathrm{\Phi }_3`$ (7)
$`i\mathrm{\Omega }_c^{}\mathrm{\Phi }_2,`$
where we have introduced the recoil velocity $`v_{\mathrm{rec}}\mathrm{}k_p/m`$. The dimensionless parameter $`\eta \mathrm{\Delta }k/k_p=1k_c^{}/k_p`$ describes the momentum transfer of the light fields to the atoms. In a degenerate $`\mathrm{\Lambda }`$ system with co-propagating probe and control field there is no momentum transfer and thus $`\eta =0`$.

Similarly we find the equation of motion of the probe field propagating parallel to the rotation in slowly-varying envelope approximation

$`\left(_t+c_zik_p\mathrm{\Omega }R\right)\mathrm{\Omega }_p(z,t)=ig\mathrm{\Phi }_2^{}\mathrm{\Phi }_1.`$ (8)
We now make use of the assumption that the control field is strong compared to the probe field and treat the interaction with the fields in lowest order of perturbation in $`\mathrm{\Omega }_p`$. In this approximation the effect of the interaction on the coupling field $`\mathrm{\Omega }_c`$ can be disregarded. Furthermore we find for the matter field $`\mathrm{\Phi }_1(z,t)=\sqrt{n}=\text{const},`$ where $`n`$ is the density of atoms. Assuming stationary conditions, the corresponding equations for $`\mathrm{\Phi }_2`$ and $`\mathrm{\Phi }_3`$ read:

$`\left[\begin{array}{cc}\mathrm{\Omega }_c& k_p\mathrm{\Omega }Ri\gamma \\ \eta k_p\mathrm{\Omega }R& \mathrm{\Omega }_c^{}\end{array}\right]\left[\begin{array}{c}\mathrm{\Phi }_3\\ \mathrm{\Phi }_2\end{array}\right]=\left[\begin{array}{c}g\sqrt{n}\mathrm{\Omega }_p\\ 0\end{array}\right]`$ (15)
$`+\left[\begin{array}{c}i\left(v_{\mathrm{rec}}\mathrm{\Omega }R\right)_z\mathrm{\Phi }_2\\ i\left(\eta v_{\mathrm{rec}}\mathrm{\Omega }R\right)_z\mathrm{\Phi }_3\end{array}\right].`$ (18)
If the complex amplitudes of the matter-wave components change sufficiently slowly, we can analytically solve eqs. (18) by treating the term on the second line as a perturbation. Substituting the corresponding result into the equation for the probe field (8) leads to the group velocity

$`v_{\mathrm{gr}}=c\mathrm{cos}^2\theta +\eta v_{\mathrm{rec}}\mathrm{sin}^2\theta .`$ (19)
One recognizes that the group velocity now contains an additional term due to the momentum transfer associated with the coherence transfer if $`\eta 0`$ Dutton-Science-2001 . In the following we will assume that $`\eta 0`$, although negative values of $`\eta `$ are possible. In the latter case the minimum group velocity attainable is zero since the medium becomes opaque as soon as $`v_{\mathrm{gr}}`$ crosses the value zero. Thus if the mixing angle exceeds the critical value $`\mathrm{tan}^2\theta _{\mathrm{crit}}c/v_{\mathrm{rec}}=mc^2/\mathrm{}\omega _p`$ the excitation propagates in the medium essentially as a matter-wave with the recoil velocity.

Under stationary conditions and keeping only terms up to first order in the rotation velocity $`\mathrm{\Omega }`$, one finds the simple probe-field equation

$`_z\mathrm{ln}\mathrm{\Omega }_p(z)=i{\displaystyle \frac{2\pi }{\lambda c}}\mathrm{\Omega }R{\displaystyle \frac{\xi (z)+\eta \frac{mc^2}{\mathrm{}\omega _p}}{\xi (z)+\eta }},`$ (20)
where $`\lambda `$ is the probe-field wavelength and we have introduced the parameter

$`\xi (z){\displaystyle \frac{\mathrm{tan}^2\theta _{\mathrm{crit}}}{\mathrm{tan}^2\theta (z)}}={\displaystyle \frac{v_{\mathrm{gr}}(z)}{v_{\mathrm{rec}}}}\eta ,`$ (21)
where in the last equation $`v_{\mathrm{gr}}c`$ was assumed. (Note that $`v_{\mathrm{gr}}\eta v_{\mathrm{rec}}`$.) Eq.(20) describes a simple phase shift without absorption losses due to perfect EIT. Two counter-propagating probe fields will thus experience the Sagnac shift

$`\mathrm{\Delta }\varphi `$ $`=`$ $`{\displaystyle \frac{2\pi \mathrm{\Omega }R}{\lambda c}}{\displaystyle dz\frac{\xi (z)}{\xi (z)+\eta }}`$ (22)
$`+{\displaystyle \frac{\mathrm{\Omega }R}{\mathrm{}/m}}{\displaystyle dz\frac{\eta }{\xi (z)+\eta }}.`$
The integration over $`z`$ takes into account that the group velocity can be different in different parts of the interferometer loop (see Fig. 1). If $`\xi \eta mc^2/\mathrm{}\omega _p`$, i.e. for $`v_{\mathrm{gr}}v_{\mathrm{rec}}mc^2/\mathrm{}\omega _p`$, eq.(22) reproduces the Sagnac phase shift of an optical gyroscope, eq.(1). On the other hand if $`\eta 0`$ and $`\xi \eta `$ (22) approaches the matter-field phase shift. Furthermore one recognizes that in the absence of momentum transfer, i.e. for $`\eta =0`$, the phase shift is for all values of $`\xi `$ identical to that of a light interferometer. Thus the reduction of the group velocity alone does not affect the Sagnac phase shift. One should note that this is in contrast to the strong Fresnel dragging always present in a medium with a small group velocity. The difference to Fresnel dragging arises because in the Sagnac interferometer there is no motion of the medium relative to the source of the waves and consequently no simple Doppler-shift. In order to distinguish rotational phase shifts from those resulting from linear accelerations it is necessary to use a symmetric interferometer setup. Since the matter wave is generated by a coherent transfer from light we can combine the rotational sensitivity per unit area with a large interferometer area utilizing well established guiding techniques for light. This is the main result of the present paper. The enhancement of the Sagnac phase shift due to an atomic vapor with EIT put into the beam path of a given optical interferometer is shown in Fig. 2. For simplicity it is assumed that the vapor fills the whole volume of the optical beam path.

In order to estimate the potential enhancement of rotational sensitivity we now have to discuss the limitations on the minimum value of $`\xi `$, which arise mainly from two-photon Doppler-shifts and velocity changing collisions. As seen above, an enhancement of the rotational sensitivity requires $`\eta >0`$. In this case there is a two-photon Doppler-shift caused by thermal motion of atoms in the rotating frame, which leads to a finite absorption of the medium. As the absorption grows with decreasing group velocity, this represents an essential limitation to the hybrid interferometer. To discuss the effects of a thermal velocity distribution we decompose the three matter fields into velocity classes according to

$`\mathrm{\Phi }_\mu (z,t)`$ $`=`$ $`{\displaystyle \underset{v}{}}\mathrm{\Phi }_\mu ^v(z,t)\mathrm{e}^{i(qz\nu t)}`$ (23)
where $`\mathrm{}q=mv`$ is the initial velocity of the atoms and $`\mathrm{}\nu =mv^2/2mv\mathrm{\Omega }R`$. We then can proceed as above and obtain the same equations of motion for the slowly-varying velocity components $`\mathrm{\Phi }_\mu ^v`$ as in eqs.(18) except for the replacement $`\mathrm{\Omega }R\mathrm{\Omega }Rv`$.

Keeping again only linear contributions of the angular velocity $`\mathrm{\Omega }`$ to the phase, we find the following propagation equation for the probe field

$`_z\mathrm{ln}\mathrm{\Omega }_p(z)=ik_p{\displaystyle \frac{\mathrm{\Omega }R}{c}}{\displaystyle \frac{\xi (z)+\eta \frac{mc^2}{\mathrm{}\omega _p}}{\xi (z)+\eta }}`$
$`2ik_p\left({\displaystyle \frac{\gamma \eta k_pc}{g^2n}}{\displaystyle \frac{\overline{v^2}}{v_{\mathrm{rec}}^2}}{\displaystyle \frac{1}{\xi (z)\left(\xi (z)+\eta \right)}}\right)^2`$ (24)
$`k_p{\displaystyle \frac{\gamma \eta k_pc}{g^2n}}{\displaystyle \frac{\overline{v^2}}{v_{\mathrm{rec}}^2}}{\displaystyle \frac{1}{\xi (z)\left(\xi (z)+\eta \right)}},`$
which has the simple solution

$`\mathrm{\Omega }_p(z)=\mathrm{\Omega }_p(0)\mathrm{e}^{i\varphi (z)}\mathrm{e}^{\kappa z}.`$ (25)
In eq. (24) $`\overline{v^2}`$ denotes the mean square of the velocity distribution with $`\overline{v}=0`$.

The first term in eq. (24) gives the Sagnac phase shift as derived previously, the second term is a constant phase contribution which however cancels in the difference phase of the two counter-propagating waves and is thus without consequence. The important new contribution due to the atomic velocity distribution is the third term in (24) which gives rise to absorption which grows with decreasing $`\xi `$. Thus $`\xi `$, resp. the group velocity, cannot be made arbitrarily small since the growing losses will lead to a decreasing signal-to-noise ratio. Whether or not the matter-wave limit $`\xi \eta `$ can be reached will depend on the velocity spread or temperature of the atoms. In deriving eq. 24 we have made furthermore the following simplifying assumptions $`\eta k_p^2\overline{v^2}|\mathrm{\Omega }_c|^2,|\mathrm{\Omega }_c|^4/\gamma ^2.`$ These are however well justified as long as the total amplitude decrease due to absorption is less than e-1 and $`\xi \alpha ^1`$, where $`\alpha g^2nz/\gamma c`$ is the opacity of the medium in the absence of EIT, which is typically in the range between $`1`$ and $`10^2`$. Fig. 3 shows the absorption coefficient

$`\kappa L={\displaystyle \frac{(k_pL)^2}{\eta \alpha }}{\displaystyle \frac{T}{T_{\mathrm{rec}}}}{\displaystyle \frac{1}{\xi (\xi +1)}}.`$ (26)
for three different temperatures $`T`$ of the atomic sample relative to the recoil temperature $`T_{\mathrm{rec}}`$. In the left figure $`\alpha =100,\lambda =500`$nm, and $`L=100\mu `$m, which corresponds to a typical situation in an ultra-cold atomic gas in a trap. In the right figure $`\alpha =10,\lambda =500`$nm, and $`L=1`$cm, corresponding to a typical vapor cell. One recognizes that by approaching the matter-wave regime the absorption experiences a steep increase at a certain value of $`\xi `$. Thus in order to retain a reasonable signal-to-noise ratio, i.e. for absorption coefficients of $`\kappa L1`$, there is a minimum value of $`\xi `$ (resp. a minimum value of $`v_{\mathrm{gr}}`$). Fig. 3 shows that it is only possible to take advantage of the maximum possible sensitivity enhancement of $`10^{11}`$ if the atomic sample is cooled down to the recoil limit or even below. Nevertheless it is still possible to gain a few orders of magnitude compare to the pure light regime by moderate cooling ($`T10^3T_{\mathrm{rec}}`$).

Another important limitation to the minimum value of $`\xi `$ is set by the decoherence of motional degrees of freedom of the matter-wave component caused by collisions. Although superpositions of internal states can survive very many collisions, in particular if they are hyperfine components of the electronic ground state, velocity kicks quickly destroy the coherence between motional states. To account for the latter in a quantitative way a kinetic theory needs to be developed. Such a theory goes however beyond the scope of the present paper and will be the subject of future work. We here can only give an upper bound for the corresponding limit on $`\xi `$ (or equivalently the group velocity $`v_{\mathrm{gr}}`$) by requiring that the pulse delay time $`\tau _\mathrm{d}`$ is shorter than the average time $`\tau _{\mathrm{coll}}`$ between two successive velocity-changing collisions:

$`\left({\displaystyle \frac{v_{\mathrm{gr}}}{v_{\mathrm{rec}}}}\right)_{\mathrm{min}}{\displaystyle \frac{L}{v_{\mathrm{rec}}\tau _{\mathrm{coll}}}}=Ln\sigma \sqrt{{\displaystyle \frac{T}{T_{\mathrm{rec}}}}}`$ (27)
where $`\sigma `$ is the collisional cross section. E.g. in a gas cell of $`L=1`$cm, $`n=10^{11}`$cm-3, assuming $`\sigma =10^{12}\mathrm{..10}^{10}`$cm2 yields $`v_{\mathrm{gr}}^{\mathrm{min}}/v_{\mathrm{rec}}0.1\mathrm{}10\sqrt{T/T_{\mathrm{rec}}}`$. When a coherent sample of atoms such as a BEC is used, the limitation due to collisions is of course absent.

In the present paper we have shown that it is possible to combine the advantages of light and mater-wave gyroscopes, to measure the Sagnac phase-shift making use of the coherence and momentum transfer associated with ultra-slow light propagation in cooled atomic vapors.

M.F. would like to acknowledge discussions with L.J. Wang, R. Chiao and U. Leonhardt during the KITP workshop on “Fast and slow light and meta-materials” at the UC Santa Barbara as well as with I. Bialynicki-Birula.