\section{ Introduction}

One of the most challenging problems in modern physics is the unification of the gravitational
interaction with other known interactions in nature. Many attempts
 involve going to higher dimensions and postulating the existence of extra
spatial dimensions.
If these  extra dimensions really exist, one must explain why they are not seen. The usual
answer is that they curl into an extremely small compactified manifold, possibly
   as small as the Planck length,  $l_{pl}=1.6\times 10^{-33}$ cm. Therefore
low-energy physics should be
insensitive to them until distances of the compactification scale are being
probed. In general, one has the possibility of observing the presence of the
extra dimensions in a scattering experiment in which energies greater than that
associated with the compactification scale are achieved. Many extra dimension
models are constained to have extra dimensions no larger than about
$ (1 {\rm Tev})^{-1} = 2 \times 10^{-17} {\rm cm}$.
 However, if only gravity
propagates in the extra dimensions, the upper bound can be much larger.  A
recent proposal is that the fundamental scale of quantum
gravity  can be as low as few Tev and the observed weakness of gravity is the
result of large extra dimensions in which only gravity can propagate
\cite{ADD,ADD2,ANT}.
  The  size of extra dimensions (of which there must be at least
two) can be as large as  1 mm in this type of model.

 However, a question arises naturally as to whether there are any lower bounds
on the sizes of extra dimensions.
 It is the common belief that the existence of extra dimensions has no effect
on low-energy physics
as long as  they are extremely small.
We will argue in this letter that this need not be the case.
The reason is lightcone fluctuations arising from the quantum
gravitational vacuum fluctuations, due to compactification of spatial
dimensions\cite{YUF,Ford95}.  The compactification of spatial
dimensions gives rise to stochastic fluctuations in the apparent speed of light
which are in principle observable.  Basically,  the smaller the size of the
compactified dimensions , the larger are the fluctuations that  result.
This is closely related to the Casimir effect,  the vacuum energy occuring
whenever boundary conditions are imposed on a quantum
field.  The gravitational Casimir energy in the five-dimensional case with one
compactified
spatial dimension was studied in \cite{APCH}, where a nonzero energy density
was found, which
tends to make the extra dimension  contract. This raises the question of
stability of the extra
dimensions. It is possible, however, that the Casimir energy arising from the
quantum
gravitational field and other matter fields may be made to cancel each
other\cite{RRT,Tsokos}, thus stabilizing the extra dimensions. Quantum lightcone
fluctuations  due to the compactfication of
 spatial dimensions\cite{YUF}, although similar in nature to the Casimir effect,
 come solely from gravitons. Hence, no similar cancelation is to be
 expected.

Note that all of these quantum effects increase in magnitude as the
compactification scale decreases. The effect of compatification
in, for example, a five dimensional model can be thought of as producing
an infinite tower of massive modes with masses inversely proportional
to the compactification scale. One might think that the contribution of these
modes to radiative correction would decrease with decreasing
compactification scale, whereas in fact the reverse is
actually the case. Detailed calculations of the vacuum polarization \cite{F80}
and of the electron self-energy \cite{Y79} have been performed which reveal
growing effects with decreasing compactification length. This can be
understood as a consequence of the uncertainty principle; the fluctuations of
quantum fields confined in a finite region increase as the size of the
region decreases.

In a recent work, we studied the light cone fluctuations in  four
dimensional flat spacetime
with compactification in one spatial dimension\cite{YUF}. It was found that
 these fluctuations,
although typically of the order of the  Planck scale, can get larger for path
lengths large compared
to the compactification scale. In particular, the mean deviation from the
classical
propagation time, $\Delta t$, is proportional to the square root of the travel
distance, $r$.
 In this paper
we apply the formalism to  spacetimes with extra dimensions, and
examine, in particular,  the five dimensional case
 to demonstrate a possible observable
consequence of  compactification of the extra dimension.
  To begin, let us examine a $d$ dimensional spacetime with $d-4$  extra
dimensions.
 Consider a flat background  spacetime  with a linearized perturbation
$h_{\mu\nu}$ propagating upon it , so the spacetime metric
may be written as  $
ds^2  = (\eta_{\mu\nu} +h_{\mu\nu})dx^\mu dx^\nu
= dt^2 -d{\bf x}^2 + h_{\mu\nu}dx^\mu dx^\nu \, ,  $
where the indices $\mu,\nu$ run through $0,1,2,3,...,d-1$.
Let $\sigma(x,x')$ be one half of the squared geodesic distance between
 a pair of spacetime points $x$ and $x'$,   and $\sigma_0(x,x')$
be the corresponding quantity in the flat background.
In the presence of a linearized metric  perturbation
 $h_{\mu\nu}$, we may expand $
\sigma = \sigma_0 + \sigma_1 + O(h^2_{\mu\nu}) \, .$
Here $\sigma_1$ is first order  in $h_{\mu\nu}$.
If we quantize  $h_{\mu\nu}$,
then  quantum gravitational vacuum
fluctuations will  lead to fluctuations in the geodesic separation, and
therefore induce
 lightcone fluctuations.  In particular, we have
$\langle \sigma_1^2 \rangle \not= 0$, since $\sigma_1$  becomes a quantum
operator when the metric perturbations are quantized.   The quantum lightcone
fluctuations give rise to
stochastic fluctuations in the speed of light, which may produce an observable
time delay or advance $\Delta t$ in the arrival times of pulses. Note that
this model uses a linearized approach to quantum gravity which is expected
to be a limit of a more exact theory. In the absence of a full theory,
this seems to be the most conservative way to compute quantum gravity effects.
One might comtemplate doing a one-loop calculation of an S-matrix element,
along the lines of those in Refs.~\cite{F80,Y79} for electrodynamics.
However, one would need to find a way to deal with the nonrenormalizability
of one-loop quantum gravity coupled to other fields \cite{DvN}.

Let us  consider the propagation of light pulses between a source and a
detector separated by a distance $r$ on a flat  background with
quantized linear perturbations. In Ref.~\cite{Ford95} it was shown that
the root-mean-squared fluctuation in the propagation time is
\begin{equation}
\Delta t = {\sqrt{\langle \sigma_1^2 \rangle_R}\over r} \, , \label{eq:delt1}
\ee
where $\langle \sigma_1^2 \rangle_R$ is a renormalized expectation value, which
was assumed to be positive. We can give an alternative derivation which applies
in the case $\langle \sigma_1^2 \rangle_R < 0$.
For a pulse which is delayed or advanced by
time $\Delta t$,
which is much less than $r$, one finds
\begin{equation}
\sigma=\sigma_0+\sigma_1+....={1\over2}[(r+\Delta t)^2-r^2]\approx r\Delta t\,.
\ee
Take the fourth power of  the above equation and average over a given  quantum
state of gravitons
$|\phi\rangle$ (e.g. the vacuum states associated with  compactification of
spatial dimensions),
\begin{equation}
\Delta t^4={\langle \phi| \sigma_1^4 |\phi\rangle\over r^4}\,.
\label{eq:TPHI}
\ee
This result is, however,  divergent due to the formal divergence of
 $\langle\phi| \sigma_1^4 |\phi\rangle $.
 We can define $\langle \phi| \sigma_1^4 |\phi\rangle$ by normal ordering,
and let
 $
\langle \phi| \sigma_1^4 |\phi\rangle_R=
\langle \phi|:\sigma_1^4:|\phi\rangle\,$.
For a free field $\psi$, Wick's theorem yields
$
\langle:\psi^4:\rangle=3\langle:\psi^2:\rangle^2=3\langle \psi^2 \rangle_R ^2
$, where the expectation value is in the vacuum state.
Hence a suitable meaure of the  deviation from the classical propagation time is
given by
\begin{equation}
\Delta t = {(3\langle \sigma_1^4 \rangle_R)^{1/4}\over r}={3^{1/4}
\sqrt{|\langle \sigma_1^2 \rangle_R|}\over r}
\approx{\sqrt{|\langle \sigma_1^2 \rangle_R|}\over r}\,.
\label{eq:MDT}
\ee
Apart from small numerical factors which we ignore, this is equivalent to
replacing $\langle \sigma_1^2 \rangle_R$ by $|\langle \sigma_1^2 \rangle_R|$
in Eq.~(\ref{eq:delt1}).
The gauge invariance of this expression has been analyzed recently\cite{YUF}.

Note, however, that $\Delta t$ is the ensemble averaged deviation, not
necessarily the expected variation in
flight time, $\delta t $,  of two pulses emitted close together in time.
The latter is given by   $\Delta t$ only when two
successive pulses are uncorrelated.
This point is discussed in detail in Ref.~\cite{Ford96}.
These stochastic fluctuations in the apparent velocity of light arising from
quantum gravitational
fluctuations are in principle observable, since
they may lead to a spread in the arrival times of pulses from
distant astrophysical sources, or the broadening of the spectral lines.
This can be used to place a lower bound on the size
of the extra dimension. Lightcone fluctuations and their possible astrophysical
observability have been recently discussed in a somewhat different framework
in Refs.~\cite{NvD,AEMN,EMN}.

  In order to find $\Delta t$ in a particular situation,
we need to calculate the quantum expectation value  $\langle \sigma_1^2
\rangle_R$ in any chosen quantum state, which can be shown to
be given by  \cite{YUF,Ford95}
\begin{equation}
\langle \sigma_1^2 \rangle_R ={1\over 8}(\Delta r)^2
\int_{r_0}^{r_1} dr \int_{r_0}^{r_1} dr'
\:\,  n^{\mu} n^{\nu} n^{\rho} n^{\sigma}
\:\, G^{(1)R}_{\mu\nu\rho\sigma}(x,x') \,.
\label{eq:interval}
\ee
Here $ dr=|d{\bf x}|$,  $\Delta r=r_1-r_0$ and $ n^{\mu} =dx^{\mu}/dr$. The
integration is taken along the null geodesic connecting two points $x$ and
$x'$, and
 $G^{(1)R}_{\mu\nu\rho\sigma}(x,x')$ is the graviton Hadamard  function,
 understood to be suitably renormalized.
We shall now work  in the transverse-tracefree
 gauge in which the gravitational perturbations have only spatial components
$h_{ij}$.
This is a gauge which retains only physical degrees of freedom.  The quantized
field operator may be expanded as
\begin{equation}
h_{ij} = \sum_{{\bf k},\lambda}\, [a_{{\bf k}, \lambda} e_{ij} ({{\bf k},
\lambda})
  f_{\bf k} + H.c. ].
\end{equation}
Here H.c. denotes the Hermitian conjugate, $\lambda$ labels the ${1\over 2}
(d^2-3d)$ independent
 polarization states,   $f_{\bf k}$  is the mode function,
 and the $e_{\mu\nu} ({{\bf k}, \lambda})$
are polarization tensors. (Units in which $32\pi G_d =1$, where $G_d$ is
Newton's constant in d dimensions  and in which $\hbar =c =1$  will be used in
this paper.)
Now suppose  the $(d-1)$-th dimension  is compactified into a
small size of
 periodicity length $L$,
so the mode function is given by
$
f_{\bf k} = (2\omega (2\pi)^{d-2}L)^{-{1\over 2}}  e^{i({\bf k \cdot x} -\omega
t)}
	\label{eq:mode2}
$
with
$k_{d-1}={2\pi n\over L}, n=0,\pm 1, \pm 2, \pm 3,...\,.$
Let us denote the associated vacuum state by $|0_L\rangle$.
 In order to calculate the
gravitational vacuum fluctuations due to compactification of the extra
 dimension, we need the
renormalized graviton Hadamard function with respect to the vacuum state
$|0_L\rangle$,
$G^{(1)R}_{ijkl}(x,x')$, which can be seen to be given by an image sum of the
corresponding Hadamard function for the uncompactified Minkowski
 vacuum, $G_{ijkl}^{(1)}$ :
\begin{equation}
 G_{ijkl}^{(1) R}(t, x_{d-1},\,t',x_{d-1}')
 ={\sum_{n=-\infty}^{+\infty}}^{\prime}G_{ijkl}^{(1) }(t,x_{d-1}, \,t',x_{d-1}'
+nL)\, ,
\ee
where the prime on the summation indicates that the $n=0$ term is excluded and
the notation
$(t, x_1,.., x_{d-2}, x_{d-1})\equiv(t, x_{d-1}) $ has been adopted.

  From now on, we will restrict ourselves to the five dimensional case with one
compactified dimension. Because the model discussed in this paper has only
one extra dimension, the results do not bear directly on the large extra
dimension models \cite{ADD,ADD2}, which require at least two extra dimensions.
Higher dimensional cases will be discussed in detail
elsewhere~\cite{YUF2}.
To see how light cone fluctuations arise in the usual uncompactified space
as a result of compactification of the extra dimension,
let us consider a light ray traveling along the
$x_1$ direction from point $a$ to point $b$, which is perpendicular to the
direction of
compactification. The relevant graviton two-point function is $G_{xxxx}$,
which can
be shown to be given by
\begin{eqnarray}
G_{xxxx}(t,x_4,\, t',x_4')&=&2[ D(t,x_4,\, t',x_4')-2F_{xx}(t,x_4,\, t',x_4')
\nonumber\\
&&+H_{xxxx}(t,x_4,\, t',x_4')]\,.
\end{eqnarray}
Here $D(x,x')$,  $F_{ij}(x,x')$ and $H_{ijkl}(x,x')$ are functions given by.
\begin{eqnarray}
 	D(x,x')&=&
{Re\over{(2\pi)^4}}\int\, {d^4{\bf k}\over{2 \omega}}e^{i{\bf k}
\cdot({\bf x}-{\bf x'})}e^{-i\omega(t-t')}\nonumber\\
&&={1\over 8\pi^2}{1\over (R^2-\Delta t^2)^{(3/2)}} \,,
\end{eqnarray}

\begin{eqnarray}
 	F_{ij}(x,x')&&={Re\over{(2\pi)^4}}\partial_i\partial'_j\int\,
{d^4{\bf k}\over{2 \omega^3}}e^{i{\bf k} \cdot({\bf x}-{\bf x'})}
e^{-i\omega(t-t')}\nonumber\\
&&={1\over{8\pi^2}}\partial_i\partial'_j\,
\left({\sqrt{R^2-\Delta t^2}\over R^2}\right)
 \,,
\end{eqnarray}
and
\begin{eqnarray}
 	H_{ijkl}(x,x')&&={Re\over{(2\pi)^4}}
   \partial_i\partial^{\prime}_j\partial_k\partial_l^{\prime}
\int\, {d^4{\bf k}\over{2 \omega^5}}e^{i{\bf k} \cdot({\bf x}-
{\bf x'})}e^{-i\omega(t-t')}
\nonumber\\
&&=0\,.
\end{eqnarray}
Here $R=|{\bf x-x'}|$ and $\Delta t=t-t'$.

 Inserting the above results into Eq.~(\ref{eq:interval}),
carrying out the differentiation in the function $F_{x_1x_1}$, using the fact
that $x_1-x_1'=\Delta t$,  and then performing the integration, we finally find
\begin{equation}
\langle \sigma_1^2 \rangle_R ={r^2\over 32\pi^2L}{\sum_{n=1}^{\infty}}\,
\left [
{8\over n}\ln\bigl(1+{\rho^2\over n^2}\bigr)-{2\rho^2\over n^3}-{8\rho^2\over
(\rho^2+n^2)n}\right]\,.
\ee
Here we have defined $r=a-b$ and a dimensionless parameter $\rho=r/L$.  We are
interested  in the case in which
$\rho \gg 1$. It then follows that the summation is dominated, to the leading
order, by
the second term,
\begin{equation}
\langle \sigma_1^2 \rangle_R \approx -
{r^2\over 16\pi^2L} {\sum_{n=1}^{\infty}}\,{\rho^2\over n^3}
 = -{\zeta(3) r^2\rho^2\over 16\pi^2L}\,,
\ee
where $\zeta(3)$ is the Riemann-zeta function.
So, the mean deviation from the classical propagation time due to the lightcone
fluctuations
is
\begin{equation}
\Delta t \approx \sqrt{ {\zeta(3)\over 16\pi^2L}}\sqrt{32\pi G_5}\,\rho
\approx \biggl({ r\over L}\biggr)\,t_{pl}\,.
\label{eq:t1}
\ee
Here we have used the fact that $G_5=G_4\, L$, and   $t_{pl} \approx 5.39\times
10^{-44}
s$ is the Planck time.

Here $\Delta t$ increases linearly with the propagation
distance, in contrast to the square root growth found in four dimensional
compactified spacetime \cite{YUF}.
Equation~(\ref{eq:t1}) also reveals that the lightcone fluctuation effect is
inversely related to $L$, the compactification length. This is due to the
increased quantum fluctuations of fields confined in a smaller region. When
$r$ is of cosmological dimensions and $L$ is sufficiently small, the effect
is potentially observable.

Before we proceed further, it should be noted that we have set
$\langle \sigma_1^2 \rangle_R =0$ when $L \rightarrow \infty$. This is the
most natural choice of renormalization, corresponding to the effect of the
graviton fluctuations vanishing in the limit of noncompactified spacetime.
This is analogous to setting a Casimir energy density to zero in the limit
of infinite plate separation.
There is, however, another logical possibility.
This is to set $\langle \sigma_1^2 \rangle_R =0$ at a finite value of $L$.
If $L$ is constant, then this procedure removes the lightcone fluctuations
and sets $\Delta t =0$. In our view, this is an unsatisfactory solution,
as there seems to be nothing in the theory which picks out a particular
finite value of $L$. In any case, changes in $L$ would still be observable,
as one could at most set $\langle \sigma_1^2 \rangle_R =0$ at one point
along the path of a light ray. Let $L_i$ be the compactification length at
the time of emission, and $L_f = L_i(1 +\delta)$, where $|\delta| \ll 1$,
be that at the time of detection. Then the variation in flight times must
satisfy~\cite{YUF2}
\begin{equation}
\Delta t \agt \sqrt{|\delta|}\, \biggl({ r\over L}\biggr)\,t_{pl} \,.
\label{eq:Tvari}
\ee
The limits on the time variation of fundamental constants place various upper
bounds \cite{KPW,JDB} on $|\delta|$ in the range between $10^{-2}$
and $10^{-10}$.

The variation in the flight time of pulses, $\Delta t$, can apply to the
successive wave crests of a plane wave. This leads to a broadening of spectral
lines from a distant source. Note, however, that  $\Delta t$ is the
 expected variation in the arrival times of two successive pulses only when
they are uncorrelated.
Following
Ref.~\cite{Ford96}, we have examined~\cite{YUF2} the correlation between
two successive pulses separated in time by $T$, and found that if $r \gg T, L$,
the
two pulses are uncorrelated when $T\gg L^2/r$. Thus the pulses can be
uncorrelated
even when $T \ll L$. Now
suppose that the
 experimental fractional resolution for a particular spectral line of period
$T$ is $\Gamma$. Then we must have ${\Delta t/ T}\leq \Gamma$, which leads
to a lower bound on $L$ of
\begin{equation}
L\agt {r\,t_{pl}\over \Gamma T} = \frac{4\times 10^4 cm}{\Gamma} \,
\left(\frac{r}{1000 Mpc}\right)\, \left(\frac{E}{1 MeV}\right) \,,
\label{eq:Lbound}
\ee
where $E$ is the energy of a photon with period $T$.
This bound can be trusted only when the condition for uncorrelated pulses,
\begin{equation}
L\ll \sqrt{rT} = 6\times 10^8 cm\, \left(\frac{r}{1000 Mpc}\right)^\frac{1}{2}
\, \left(\frac{1 MeV}{E}\right)^\frac{1}{2}
\,,  \label{eq:uncorr}
\ee
is satisfied.

The strongest lower  bound would  be deduced from the highest frequency
spectral  lines of the most distant sources, with the highest observed
resolution. The best compromise between these various requirements seems to be
data from gamma ray burst sources, which involve cosmological distances and
high frequencies, albeit low resolution. The use of gamma ray burst to
constrain quantum gravity models was discussed by Amelino-Camelia, {\it et al}
\cite{AC}. A typical burst is GRB990123, which involved\cite{GRB990123}
gamma rays at an energy of the order of $1\, MeV$ from a source
with a redshift of at least $z = 1.6$. Assuming a Hubble constant of
$H_0 = 65 km/s/Mpc$ and a matter dominated universe, this corresponds to a
distance of $r \agt 2400\, Mpc$. If we take the resoluton $\Gamma$ to be of
order unity, this leads to
a lower bound on $L$ from Eq.~(\ref{eq:t1}) of $L \geq 10^5 cm$. This lower
bound satisfies Eq.~(\ref{eq:uncorr}). This is a remarkably strong bound
which would seem to rule out the five dimensional theory. Even if one adopts
the approach of only using Eq.~(\ref{eq:Tvari}) to bound $|\delta|$, the
result is $|\delta| \leq 10^{-10} (L/ 1\, cm)^2$, which for $L \ll 1 cm$
is much stronger than
the bounds cited above based upon time dependence of fundamental constants.
It is of interest to note that reasonably strong bounds on $L$ can be
obtained from much lower frequency sources. Microwave lines from quasars
\cite{YUF2}, and the cosmic microwave background \cite{DFY} both yield lower
bounds on $L$ of the order of a few tenths of a millimeter.

To conclude, we have demonstrated, in the case of one extra dimension, that
the large quantum lightcone fluctuations due to the compactification of the
extra dimension require the size of the extra dimension to be macroscopically
large. This result seems to rule out the five dimensional Kaluza-Klein theory,
or at the very least, place strong limits on the rate of change of the extra
dimension. We must point
out that the rate of growth of $\Delta t$ with $r$  depends crucially on the
number and nature of the
spatial dimensions. In four dimensions, $\Delta t\propto \sqrt{r/L}$,
while in five dimensions
$\Delta t \propto r/L$. We have also analyzed six and higher dimensional models
in which the extra dimensions are flat. Here we find \cite{YUF2}
that $\Delta t$ grows
only as a  logarithmic function of $r/L$. In these models, data on spectral
lines from distant sources yield no significant constraints on the
compactification scale.

In principle, the effect discussed in this paper should be a generic
phenomenon to be expected whenever there are small extra dimensions. This
follows from the uncertainty principle argument given above, to the effect
that confining the graviton modes in a small region of space should give rise
to large fluctuations inversely related to the size of the region. However,
there is a possibility of subtle cancellations making the net effect smaller
than would naively be expected. This is what happens in the models with two or
more flat extra dimensions.
However, with two or more extra dimensions, it would
also be possible to have the extra dimensions curved. The lightcone fluctuations
in such models have not yet been examined. It is entirely possible that they
might exhibit a rate of growth intermediate between linear and logarithmic
functions. The phenomenon discussed in this
paper can not only constrain models with extra dimensions, but could
conceivably lead to positive confirmation of the existence of such dimensions.

Acknowledgment:
We would like to thank R. Di Stefano and W. Waller for useful discussions.
This work was
supported in part by the National Science Foundation under Grant PHY-9800965.


\section{Introduction} \label{sec:intro}

        QCD admits a useful and elegant expansion in powers of $1/N_c$,
where $N_c$ is the number of colors~\cite{tHooft}.  This expansion has
been utilized successfully in baryon effective field theories to isolate
the leading and subleading contributions to a variety of physical
observables~\cite{revs}.

Here we study the mass spectrum of the nonstrange, $\ell=1$ baryons
(associated with the SU(6) {\bf 70}-plet for $N_c=3$) in a large-$N_c$
effective theory~\cite{CCGL1,CCGL2}.   We describe the states as a
symmetrized ``core'' of $(N_c-1)$ quarks in the ground state plus one
excited quark in a relative $P$ state.  ``Quarks'' in the effective theory
refer to eigenstates of the spin-flavor-orbit group, SU(6) $\times$ O(3),
such that an appropriately symmetrized collection of $N_c$ of them have
the quantum numbers of the physical baryons.  Baryon wave functions
are antisymmetric in color and symmetric in the spin-flavor-orbit indices
of the quark fields.  While this construction assures that we obtain
states with the correct total quantum numbers, we do not assume that
SU(6) is an approximate symmetry of the effective Lagrangian.  Rather,
we parameterize the most general way in which spin and flavor symmetries
are broken by introducing a complete set of quark operators that act on
the baryon states.  Matrix elements of these operators are hierarchical
in $1/N_c$, so that predictivity can be obtained without recourse to
ad hoc phenomenological assumptions.

The nonstrange {\bf 70}-plet states which we consider in this analysis
consist of two isospin-$3/2$ states, $\Delta_{1/2}$ and $\Delta_{3/2}$,
and five isospin-$1/2$ states, $N_{1/2}$, $N^\prime_{1/2}$, $N_{3/2}$,
$N^\prime_{3/2}$, and $N^\prime_{5/2}$.  The subscript indicates total
baryon spin; unprimed states states have quark spin $1/2$ and primed
states have quark spin $3/2$.  These quantum numbers imply
that two mixing angles, $\theta_{N1}$ and $\theta_{N3}$, are necessary
to specify the total angular momentum $1/2$ and $3/2$ nucleon mass
eigenstates, respectively.  Definitions for these angles can be found
in Ref.~\cite{CGKM}.

\section{Operators Analysis} \label{sec:melements}

To parameterize the complete breaking of SU(6)$\times$O(3), it is
natural to write all possible mass operators in terms of the generators
of this group.   The generators of orbital angular momentum are denoted
by $\ell^i$, while $S^i$, $T^a$, and  $G^{ia}$ represent the spin,
flavor, and spin-flavor generators of SU(6), respectively. The generators
$S_c^i$, $T_c^a$, $G_c^{ia}$ refer to those acting upon the $N_c-1$ core
quarks, while separate SU(6) generators $s^i$, $t^a$, and $g^{ia}$
act only on the single excited quark.  Factors of $N_c$ originate
either as coefficients of operators in the Hamiltonian, or through
matrix elements of those operators.  An $n$-body operator,
which acts on $n$ quarks in a baryon state, has a coefficient of
order $1/N_c^{n-1}$, reflecting the minimum number of
gluon exchanges required to generate the operator in QCD.
Compensating factors of $N_c$ arise in matrix elements if sums over
quark lines are coherent. For example, the unit operator $1$ contributes
at $O(N_c^1)$, since each core quark contributes equally in the matrix
element. The core spin of the baryon $S_c^2$ contributes to the masses
at $O(1/N_c)$, because the matrix elements of $S_c^i$ are of $O(N_c^0)$
for baryons that have spins of order unity as $N_c \to \infty$.  Similarly,
matrix elements of $T_c^a$ are $O(N_c^0)$ in the two-flavor case since
the baryons considered have isospin of $O(N_c^0)$, but the operator
$G_c^{ia}$ has matrix elements on this subset of states of $O(N_c^1)$.
This means that the contributions of the $O(N_c)$ quarks add incoherently
in matrix elements of the operator $S_c^i$ or $T_c^a$ but coherently
for $G_c^{ia}$. Thus, the full large $N_c$ counting of the matrix element
is $O(N_c^{1-n+m})$, where $m$ is the number of coherent core quark
generators. A complete operator basis for the nonstrange {\bf 70}-plet masses
is shown in Table~1\footnote{Some of these operators were considered
previously in Ref.~\cite{goity}.}. Index contractions are left implicit
wherever they are unambiguous, and the $c_i$ are operator coefficients.  The
tensor $\ell^{(2)}_{ij}$ represents the rank two tensor combination
of $\ell^i$ and $\ell^j$ given by
$\ell^{(2)}_{ij} = \frac{1}{2} \{ \ell_i , \ell_j \} - \frac{\ell^2}{3}
\delta_{ij}$.
\begin{table}[htb]
\caption{A complete operator basis, ${\cal O}_i$, $i=1\ldots 9$, for the
nonstrange {\bf 70}-plet masses.}
\label{matel}
\newcommand{\hphantom{$-$}}{\hphantom{$-$}}
\newcommand{\cc}[1]{\multicolumn{1}{c}{#1}}
\renewcommand{\tabcolsep}{2pc}
\renewcommand{\arraystretch}{1.2}
\begin{tabular}{@{}lllll}
\hline
$c_1${\bf 1} & $c_2\ell s$  & $c_3\frac{1}{N_c} \ell^{(2)} gG_c$ & & \\
$c_4(\ell s + \frac{4}{N_c+1} \ell tG_c)$ & $c_5\frac{1}{N_c} \ell S_c$ &
$c_6\frac{1}{N_c} S_c^2 $ & & \\
$c_7\frac{1}{N_c} tT_c $ & $c_8\frac{1}{N_c}\ell^{(2)} sS_c $
& $c_9\frac{1}{N_c^2} \ell^i g^{ja} \{ S_c^j,G_c^{ia} \}$ & & \\
\hline
\end{tabular}\\[2pt]
Operators $1$, $2$--$3$, and
$4$--$9$ have matrix elements  of order $N^1_c$,
$N^0_c$, and $N^{-1}_c$, respectively.
\end{table}
\section{Results} \label{sec:results}

        Since the operator basis in Table~1 completely spans
the $9$-dimensional space of observables, we can solve for the $c_i$
given the experimental data.  For each baryon mass, we assume that the
central value corresponds to the midpoint of the mass range quoted in
the {\it Review of Particle Properties}~\cite{RPP}; we take the one
standard deviation error as half of the stated range.  To determine
the off-diagonal mass matrix elements, we use the mixing angles
extracted from the analysis of strong decays given in
Ref.~\cite{CGKM}, $\theta_{N1}=0.61\pm 0.09$ and $\theta_{N3}= 3.04
\pm 0.15$.  These values are consistent with those obtained in
\cite{CC} from radiative decays.  Solving for the operator
coefficients, we obtain the values shown in Table~\ref{csolve}.

\begin{table}[htb]
\caption{Operator coefficients in GeV, assuming the complete set of
Table~\protect\ref{matel}.}
\label{csolve}
\newcommand{\hphantom{$-$}}{\hphantom{$-$}}
\newcommand{\cc}[1]{\multicolumn{1}{c}{#1}}
\renewcommand{\tabcolsep}{0.64pc}
\renewcommand{\arraystretch}{1.2}
\begin{tabular}{@{}lllllllll} \hline
$c_1$ & $c_2$ & $c_3$ & $c_4$ & $c_5$ & $c_6$ & $c_7$ & $c_8$ &
$c_9$ \\ \hline
$+0.470$ & $-0.036$ & $+0.369$ & $+0.089$ & $+0.087$ & $+0.418$ &
$+0.040$ & $+0.048$ & $+0.012$ \\
$\pm 0.017$ & $\pm 0.041$ & $\pm 0.208$ & $\pm 0.203$ & $\pm 0.157$
& $\pm 0.085$ & $\pm 0.074$ & $\pm 0.172$ & $\pm 0.673$ \\\hline
\end{tabular}\\[2pt]
\end{table}

        Naively, one expects the $c_i$ to be of comparable size. Using
the value of $c_1$ as a point of comparison, it is clear that there
are no operators with anomalously large coefficients.  Thus, we find
no conflict with the naive $1/N_c$ power counting rules.  However,
only three operators of the nine, ${\cal O}_1$, ${\cal O}_3$,
and ${\cal O}_6$, have coefficients that are statistically distinguishable
from zero!   A fit including those three operators alone is shown in
Table~\ref{3param}, and has a $\chi^2$ per degree of freedom is $1.87$.
Fits involving other operator combinations are studied in
Refs.~\cite{CCGL1,CCGL2}. Clearly, large $N_c$ power counting is not
sufficient by itself to explain the $\ell=1$ baryon masses---the
underlying dynamics plays a crucial role.
\begin{table}[htb]
\caption{Three parameter fit using operators ${\cal O}_1$, ${\cal
O}_3$, and ${\cal O}_6$, giving $\chi^2/{\rm d.o.f.}=$
$11.19/6=1.87$. Masses are given in MeV, angles in radians.}
\label{3param}
\newcommand{\hphantom{$-$}}{\hphantom{$-$}}
\newcommand{\cc}[1]{\multicolumn{1}{c}{#1}}
\renewcommand{\tabcolsep}{1.3pc}
\renewcommand{\arraystretch}{1.2}
\begin{tabular}{@{}lllllll}  \hline
               & Fit       & Exp.\       && & Fit       & Exp.\
\\ \hline
$\Delta(1700)$ & $1683$ & $1720\pm 50$ &&$N(1520)$ & $1530$ & $1523\pm
8$ \\
$\Delta(1620)$ & $1683$ & $1645\pm 30$ &&$N(1535)$ & $1503$ & $1538\pm
18$\\
$ N(1675) $ & $1673$ & $1678\pm 8$ && $\theta_{N1}$ & $0.45$ &
$0.61\pm 0.09$ \\
$ N(1700) $ & $1725$ & $1700\pm 50$ && $\theta_{N3}$ &$3.04$& $3.04\pm
0.15$ \\
$N(1650)$   &     $1663$ & $1660\pm 20$ &&  & & \\\hline
\end{tabular}\\[2pt]
Parameters (GeV): $c_1=0.461\pm 0.005$,
$c_3=0.360\pm 0.059$, $c_6=0.453\pm 0.030$
\end{table}

\section{Interpretation and Conclusions} \label{sec:concl}

We will now show that the preference in Table~2 for two nontrivial
operators, $\frac{1}{N_c} \ell^{(2)} g \,G_c$ and $\frac{1}{N_c} S_c^2$,
can be  understood in a constituent quark model with a single pseudoscalar
meson exchange, up to corrections of order $1/N_c^2$.  The argument goes
as follows:

The pion couples to the quark axial-vector current so that the
$\overline{q}q\pi$ coupling introduces the spin-flavor structure
$\sigma^i \tau^a$ on a given quark line. In addition, pion exchange
respects the large $N_c$ counting rules given in Section~2.  A single
pion exchange between the excited quark and a core quark is mapped
to the operators $g^{ia} G_c^{ja} \ell^{(2)}_{ij}$ and
$g^{ia} G_c^{ia}$, while pion exchange between two core quarks yields
$G^{ia}_c G^{ia}_c$. These exhaust the possible two-body operators that
have the desired spin-flavor structure. The first operator is one of the
 two in our preferred set. The third operator may be rewritten
\begin{equation}
2 G^{ia}_c G^{ia}_c = C_2 \cdot 1 - \frac{1}{2}T^a_c T^a_c
- \frac{1}{2} S_c^2
\label{eq:cas}
\end{equation}
where $C_2$ is the SU(4) quadratic Casimir for the totally symmetric core
representation (the {\bf 10}  of SU(4) for $N_c=3$).  Since the core
wavefunction involves two spin and two flavor degrees of freedom, and is
totally symmetric, it is straightforward to show that  $T_c^2=S_c^2$.  Then
Eq.~(\ref{eq:cas}) implies that one may exchange $G^{ia}_c G^{ia}_c$ in
favor of the identity operator and $S_c^2$, the second of the two operators
suggested by our fits.

The remaining operator, $g^{ia} G_c^{ia}$, is peculiar in that
its matrix element between two nonstrange, mixed symmetry states is given
by~\cite{CCGL1}
\begin{equation}
\frac{1}{N_c} \langle gG \rangle = - \frac{N_c+1}{16 N_c} +
\delta_{S,I}\frac{I(I+1)}{2N_c^2}  \,\,\, ,
\end{equation}
which differs from the identity only at order $1/N_c^2$.  Thus to order
$1/N_c$, one may make the
replacements
\begin{equation}
\{1 \mbox{ , } g^{ia} G_c^{ja} \ell^{(2)}_{ij}  \mbox{ , }
g^{ia} G_c^{ia} \mbox{ , } G^{ia}_c G^{ia}_c \} \Rightarrow
\{1  \mbox{ , }  g^{ia} G_c^{ja} \ell^{(2)}_{ij} \mbox{ , }
S_c^2\} \,\,\, .
\end{equation}
We conclude that the operator set suggested by the data may
be understood in terms of single pion exchange between quark
lines.  This is consistent with the interpretation of
the mass spectrum advocated by Glozman and Riska~\cite{gloz}.
Other simple models, such as single gluon exchange, do not
directly select the operators suggested by our analysis and
may require others that are disfavored by the data.


\section{Introduction}
Vortices and other topological field configurations play a fundamental
role in determining the properties of many fascinating materials and
control the physical mechanisms underlying several
applications. Examples include superfluids \cite{sonin},
superconductors \cite{brandt,farrell}, periodic solids \cite{mermin},
liquid crystals \cite{ss}, two-dimensional magnets \cite{msb},
propagating coherent light beams \cite{rss}, and perhaps even the hot
plasma that filled the very early Universe \cite{vs}. Common
aspects of the phenomenology of these diverse applications stem from
the mathematical similarity of the underlying field theories and their
classical static solutions.

Models of vortices as Brownian point particles, characterized by a
mass, mutual interactions, and damping, are commonly used to study the
behavior of superconductors and superfluids \cite{bmrefs}.  In this
picture the equation of motion for a single overdamped vortex, not
subject to any external forces, is extremely simple:
\begin{eqnarray}
\gamma \dot{\bf r}(t)&=&{\bf f}(t); \nonumber \\
\langle f_\mu(t)f_\nu(t') \rangle&=&\frac{2\gamma
k_BT}{M}\delta_{\mu,\nu}\delta(t-t'), \label{vbm}
\end{eqnarray}
where $\gamma$ is the damping coefficient, $M$ is the vortex mass,
${\bf f}(t)$ is a Gaussian thermal noise, $\mu,\nu\in\{x,y\}$, and
$\langle\cdots\rangle$ denotes canonical ensemble averaging.  It follows that
the vortex velocity distribution is also Gaussian.

This effective picture is extremely appealing---primarily because of
the drastic reduction of the number of degrees of freedom. However, to
our knowledge, it has not been shown to be derivable from the dynamics
of an underlying field theory. In this paper we study this question by
solving both analytically and numerically for the motion of the vortex
as an effective degree of freedom arising in a stochastic nonlinear
Schr\"odinger equation. It is well-known that the conventional
description of Brownian motion in terms of equations such as
(\ref{vbm}), that of a heavy particle interacting with light
scatterers, ignores the presence of hydrodynamic fluctuations. A more
complete description would allow us not only to test the validity of
the Brownian motion model for the vortex but also to compute
corrections to it and, possibly, to find new physical effects.

Our analytic approach rests on the use of a singular perturbation
expansion around a rigid vortex utilizing a low-temperature or
weak-noise expansion. We are able to derive a Fokker-Planck equation
for the single vortex distribution function which corresponds to that
of a passive scalar advected to a background flow [not just a simple
diffusion equation as would be predicted by Eqs.\ (\ref{vbm})]. In our
numerical work we are able to extract the diffusion constant for the
vortex motion which turns out to be in good agreement with the
theoretical prediction. Consistent with our theoretical analysis, the
vortex effective mass diverges logarithmically with the system
size. We also investigate the statistics of vortex flights and
demonstrate that the probability distribution function (PDF) has an
exponential tail implying nontrivial correlations in the background
thermal flow field.

In Section I we present the stochastic nonlinear Schr\"odinger
equation and discuss its physical relevance. In Section II we use the
low-temperature perturbative expansion about the vortex collective
coordinates to derive its equation of motion. Section III describes
our numerical techniques and the associated results. We end with
Section IV with a further discussion of our conclusions.

\section{The Stochastic Nonlinear Schr\"odinger Equation}

We study the motion of an isolated vortex by considering it as a
singularity of a classical stochastic nonlinear Schr\"odinger field
$\psi({\bf r},t)$ in two spatial dimensions. The equation of motion of this
model theory is $(a<0)$:
\begin{eqnarray}
(i\hbar-\gamma)\frac{\partial \psi}{\partial
t}&=&-\frac{\hbar^2}{2m}\nabla^2\psi
+a\psi+b|\psi|^2\psi+\eta({\bf r},t); \nonumber\\
\langle\eta({\bf r},t)\eta^*({\bf
r}',t')\rangle&=&2k_BT\gamma\delta(t-t')\delta({\bf r}
-{\bf r}'), \nonumber\\
\langle\eta({\bf r},t)\eta({\bf r}',t')\rangle&=&0.
\label{snls}
\end{eqnarray}
The damping coefficient $\gamma$ and the stochastic force $\eta({\bf r},t)$
model the coupling of the condensate to a heat bath, satisfying the
fluctuation-dissipation relation. Coefficients $a$, $b$ and $m$ are
assumed constant and independent of temperature; thermal effects are
fully described by the stochastic force.

The above equations (\ref{snls}) describe the phenomenology of
physically important systems such as atomic Bose condensates
\cite{stoof} and superfluid helium II \cite{p,gross}.  The same
equation (at $T=0$), amended as appropriate by the vector potential,
describes the order-parameter dynamics of type-II
superconductors \cite{schmid,gorkov}.  Although the above is known to
be strictly true only for a class of type-II superconductors with
magnetic impurities \cite{eliashberg}, it is often assumed that
Eqs.\ (\ref{snls}) apply to all type-II superconductors, unless
indicated otherwise.  We will henceforth study these equations as an
important paradigmatic field theory representing a broad range of
related physical phenomena.

It is convenient to make all quantities dimensionless by standard
substitutions.  We set $\psi = \tilde{\psi}\sqrt{-a/b}$, $\gamma =
\tilde{\gamma}\hbar$, $t=\tilde{t}\gamma/(-a)$, ${\bf r}=\tilde{\bf
r}\hbar /\sqrt{-2am}$, $\eta = \tilde{\eta}(-a)\sqrt{-a/b}$, and $k_B
T = -a\hbar^2/(2\tilde{\beta}mb)$.  Note that ${\bf r}$ is normalized
to a unit of length $\hbar/\sqrt{-2am}$, which is the Ginzburg-Landau
coherence length, usually denoted $\xi$.  With these substitutions,
and after dropping the tildes, Eqs.\ (\ref{snls}) become
\begin{eqnarray}
(i/\gamma-1)\frac{\partial}{\partial t} \psi& =& -\nabla^2 \psi
+[|\psi|^2-1]\psi+\epsilon\eta({\bf r},t); \nonumber\\
\langle \eta({\bf r},t) \eta^*({\bf r}',t')\rangle &=& \frac{2}{\beta}
\delta(t-t') \delta({\bf r}-{\bf r}'), \nonumber\\
\langle \eta({\bf r},t) \eta({\bf r}',t') \rangle &=& 0.
\label{dsnls}
\end{eqnarray}

The position of an isolated vortex, ${\bf R}(t)$, is defined via a contour
integral
\begin{equation}
\oint [\nabla\mbox{arg}\,\psi({\bf r},t)] \cdot d\,{\bf l}, \label{vortpos}
\end{equation}
which equals $2\pi$ when the integration path encloses the point
${\bf R}(t)$, and vanishes otherwise.

The Brownian motion of vortices as point particles has been argued to
hold primarily in the overdamped limit ($\gamma \rightarrow \infty$),
which is directly relevant to the superconducting case; this is also
the limit studied in this paper.  We envision a situation in which a
single vortex exists in the ground state of the field
$\psi({\bf r},t)$.
In a thin superconducting film one can in
principle achieve the same effect by placing the sample in a magnetic
field of one flux quantum per sample area.  In the latter case there
are many other factors of practical importance, such as geometrical
and point pinning, which however will not be dealt with here.  The
vortex is assumed to be located near the center of the spatial extent
of the condensate, so that boundary effects can be neglected.  The
radius of the sample $\Lambda_c$ provides a natural infrared cutoff.
At sufficiently low temperature, thermally induced vortex-antivortex
pairs will be exponentially suppressed and we may assume that there is
only one vortex in a spatially bounded superfluid film at all
times. This assumption has been verified numerically for a range of
temperatures and sample sizes studied.

\section{Perturbative Analysis of Vortex Transport}

In this Section we utilize the singular perturbation expansion of Kaup
\cite{kaup} to extract the vortex position from the stochastic
equation of motion (\ref{dsnls}).  We first expand the field as
\begin{eqnarray}
\psi({\bf r},t) &=& \psi_0({\bf X})+\epsilon\phi({\bf r},t); \label{psi}\\
\phi({\bf r},t) &=&\phi_1({\bf X},t) +\epsilon \phi_2({\bf X},t)+ \cdots \label{phi},
\end{eqnarray}
where ${\bf X}\equiv (X,Y)$ is the comoving coordinate, measured relative
to a moving reference point; ${\bf X}(t) = {\bf r}-{\bf R}^0(t)$.  The absence of
explicit time dependence in $\psi_0({\bf X})$ in the comoving reference
frame indicates that it is the static vortex field in the absence of
fluctuations.  Tomboulis \cite{et} has shown that introduction of a
collective coordinate such as ${\bf R}^0(t)$, while conserving the total
number of dynamical degrees of freedom, is a canonical transformation.
The fluctuating field $\phi({\bf r},t)$ can be thought of as a
superposition of harmonic modes in the background of the rigid vortex
$\psi_0({\bf X})$. The phonons are gapless \cite{pit}, therefore they will
be an important consideration at any finite temperature.  The variable
$\epsilon \ll 1$ plays the role of a small parameter as well as that of a
bookkeeping device to control the perturbation series \cite{gardiner}.
Another useful way to think about $\epsilon$ emerges when one absorbs $\epsilon$
into the definition of $\eta({\bf r},t)$ in Eq.\ (\ref{dsnls}): $\epsilon$ then
reappears as the square root of temperature in the
fluctuation-dissipation relation in Eqs.\ (\ref{dsnls}). Therefore,
the small-$\epsilon$ expansion is actually the {\it small temperature or
weak noise expansion.}  The equation of motion (\ref{dsnls}) is
required to hold at every order in $\epsilon$, which results in a hierarchy
of equations for $\psi_0({\bf X})$, $\phi_1({\bf X},t)$, etc.

Because of the added fluctuations in Eq.\ (\ref{psi}), the actual
vortex position ${\bf R}(t)$ as the singularity of the full dynamical field
$\psi({\bf r},t)$ does not in general coincide with the coordinate
${\bf R}^0(t)$, which by definition is the singularity of the rigidly
moving static vortex field $\psi_0({\bf X})$.  The two are related via
\begin{equation}
{\bf R}(t)= \vec{\rho}\,[{\bf R}^0(t),\epsilon\phi] ={\bf R}^0(t)+ \epsilon\int d^2{\bf r} \left.
\frac{\delta \vec{\rho}\,[{\bf R}^0(t),\phi]}{\delta \phi({\bf r})}
\right|_{\phi=0}\phi({\bf r},t)+\mbox{c.c.}+\cdots, \label{coord}
\end{equation}
where $\vec{\rho}\,[{\bf R}^0(t),\epsilon\phi]$ is a vector functional of the
fluctuation field $\epsilon\phi({\bf r},t)$ and of the reference point ${\bf R}^0(t)$.
We define
\begin{eqnarray}
\dot{{\bf R}}^0(t)&=&0 + \epsilon{\bf v}^0_1(t) + \epsilon^2{\bf v}^0_2(t) +\cdots,
\label{rdot1}
\\
\dot{{\bf R}}(t)&=&0 + \epsilon{\bf v}_1(t) + \epsilon^2{\bf v}_2(t) +\cdots,
\label{rdot2}
\end{eqnarray}
which makes explicit the fact that there is no vortex motion in the
absence of thermal fluctuations. Differentiating Eq.\ (\ref{coord})
with respect to $t$, we find at order $\epsilon$,
\begin{equation}
{\bf v}_1(t) = {\bf v}^0_1(t) + {\bf V}[{\bf R}(t),t], \label{v1}
\end{equation}
where ${\bf V}[{\bf R}(t),t] \equiv \frac{d}{dt}\int
d^2{\bf r}\left. \frac{\delta\vec{\rho}\,[{\bf R}(t),\phi]}{\delta\phi({\bf r})}
\right|_{\phi=0}\phi({\bf r},t)+\mbox{c.c.}$ is a time-dependent velocity
field. In the last expression we replaced ${\bf R}^0(t)$ with ${\bf R}(t)$, as
the two are equal to ${\cal O}(\epsilon)$. The statistical properties of ${\bf
V}[{\bf R}(t),t]$ are given by those of the excitation field
$\phi({\bf r},t)$\/, i.e., the phonons, and will therefore depend on
temperature.

Substituting Eqs.\ (\ref{psi}), (\ref{phi}), (\ref{rdot1}) and
(\ref{rdot2}) into Eqs.\ (\ref{dsnls}) and collecting powers of $\epsilon$
we obtain, at order $\epsilon^0$,
\begin{equation}
0=-\nabla_{\bf X}^2\psi_0-[1-|\psi_0|^2]\psi_0. \label{lowest}
\end{equation}
The properties of the static vortex solution are well known;
$\psi_0({\bf X}) = f(|{\bf X}|)\exp[i\,\mbox{arg}(X+iY)]$, where $f(x)\sim x$
for $x\ll 1$ and $f(x) \sim 1$ for $x\gg 1$. At order $\epsilon$ we find
\begin{equation}
\frac{\partial}{\partial t}\left( \begin{array}{c}
\phi_1\\ \phi^*_1 \end{array} \right) =
{\cal L} \left( \begin{array}{c} \phi_1\\ \phi^*_1 \end{array}
\right)+{\bf v}^0_1\cdot\nabla_{\bf X}\left(
\begin{array}{c} \psi_0\\ \psi^*_0 \end{array} \right) -
\left( \begin{array}{c} \eta\\ \eta^* \end{array} \right),
\label{next}
\end{equation}
where ${\cal L}$ is a Hermitean matrix
\begin{equation}
{\cal L} \equiv
\left( \begin{array}{lr}
\nabla_{\bf X}^2+[1-2|\psi_0|^2] & \psi_0^2\\
(\psi_0^*)^2 & \nabla_{\bf X}^2+[1-2|\psi_0|^2]
\end{array}
\right).
\end{equation}
The noise $\eta({\bf X},t)$ remains white in space and time also in the
moving reference frame, just as in the original Eqs.\ (\ref{dsnls}).
Differentiating Eq.\ (\ref{lowest}) with respect to $X, Y$ we find
\begin{equation}
{\cal L}\left( \begin{array}{c}
\partial_X \psi_0\\
\partial_X \psi^*_0
\end{array} \right)={\cal L}\left(
\begin{array}{c}
\partial_Y \psi_0\\
\partial_Y \psi^*_0
\end{array} \right)=0.
\end{equation}
Hence the null space of the linear operator ${\cal L}$ is spanned by the
eigenvectors
$(\partial_X \psi_0,\partial_X \psi^*_0)^T$ and
$(\partial_Y \psi_0,\partial_Y \psi^*_0)^T$.
These generate uniform translations of $\psi_0$ in the $xy$ plane and
correspond to the Goldstone modes.

It turns out \cite{kaup} that all secular terms in the
$\epsilon$\/-expansion vanish when we demand that the arbitrary perturbation
$(\phi,\phi^*)^T$ be orthogonal to the null space of the operator
${\cal L}$.  With this condition in place [which also uniquely specifies the
reference point ${\bf R}^0(t)$] we now take the scalar product of
Eq.\ (\ref{next}) with the eigenvector $(\partial_X \phi_0, \partial_X
\phi^*_0)^T$, and integrate over the sample area to obtain
\begin{equation}
0=\mbox{Re} \left\{ \int d^2{\bf X}\, (\partial_X\psi^*_0)[{\bf v}^0_1(t)
\cdot\nabla_{\bf X}\psi_0-\eta({\bf X},t)]
\right\}. \label{next2}
\end{equation}

The same equation is obtained if we use the other eigenvector in the
scalar product, only now $\partial_X\psi^*_0$ is replaced by
$\partial_Y\psi^*_0$. Solving Eq.\ (\ref{next2}) for ${\bf v}^0_1(t)$
leads to
\begin{equation}
{\bf v}^0_1(t)={\bf F}(t), \label{stoch}
\end{equation}
where the stochastic force ${\bf F}(t)$ is Gaussian with zero mean and
autocorrelator
\begin{equation}
\langle F_\mu(t)F_\nu(t')\rangle = \frac{2}
{\beta M(\Lambda_c)}
\delta(t-t')\delta_{\mu\nu},
\end{equation}
where $\mu,\nu \in \{X, Y\}$, and
\begin{equation}
M(\Lambda_c) \equiv \int^{\Lambda_c}d^2{\bf X}\,|\nabla_{\bf X} \psi_0({\bf X})|^2
\simeq 2\pi \ln\Lambda_c \label{mass}
\end{equation}
can be interpreted as the (cutoff-dependent) {\it inertial mass} of
the vortex. The fully dimensional form of this equation is
$M(\Lambda_c)\simeq (2\pi\gamma^2/b)\ln(\Lambda_c/\xi)$.  This
quantity, in a different derivation, has received the same
interpretation by \v{S}im\'{a}nek \cite{simanek}.
The so-called vortex core mass\cite{suhl}
is implicitly present in the first part of Eq.\ (18) and makes a constant
contribution to $M(\Lambda_c)$, whereas the asymptotic expression
in the latter part of Eq.\ (18) contains the
$\Lambda_c$\/-dependence which is dominant in the limit of large
$\Lambda_c$.
In charged superfluids there is yet another contribution to the vortex
inertial mass, the so-called electromagnetic mass\cite{suhl}, which
originates from the energy of the electric field generated by a moving
vortex. This contribution is absent in Eq.\ (18), which strictly applies
only to neutral superfluids.

Substituting for ${\bf v}^0_1(t)$ in Eq.\ (\ref{stoch}) from
Eq.\ (\ref{v1}) and using definitions (\ref{rdot1}) and (\ref{rdot2}),
we arrive at an equation of motion for the vortex, valid to ${\cal O}(\epsilon)$,
\begin{equation}
\frac{d}{dt}{\bf R}(t) =\epsilon{\bf F}(t)+\epsilon{\bf V}[{\bf R}(t),t].
\end{equation}
As the last step we set $\epsilon=1$ in this equation.

If we assume that ${\bf V}[{\bf R}(t),t]$ is slowly varying relative to the
noise ${\bf F}(t)$, then we may treat ${\bf V}[{\bf R}(t),t]$ as an
external field \cite{parisi} and arrive immediately at the
Fokker-Planck equation \cite{gardiner} for the density $c({\bf r},t) \equiv
\langle \delta[{\bf r}-{\bf R}(t)]\rangle$:
\begin{equation}
\frac{\partial c({\bf r},t)}{\partial t} +\nabla\cdot[{\bf V}({\bf r},t)
c({\bf r},t)] = \sigma \nabla^2 c({\bf r},t), \label{fpe}
\end{equation}
where the diffusion coefficient
\begin{equation}
\sigma \equiv \frac{1}{\beta M(\Lambda_c)}. \label{diffcon}
\end{equation}
This equation is of the same form as that for a {\it passive scalar
advected to background fluid flow} ${\bf V}({\bf r},t)$ \cite{bj,pre}. It
is well established---particularly in the best-studied case of
incompressible fluid flow $\nabla\cdot{\bf V}({\bf r},t)=0$---that in
contrast to pure diffusion, an advected passive scalar may display
{\it non-Gaussian} (for example, exponential) tails in the probability
density function (PDF) $c({\bf r},t)$. These tails arise due to nontrivial
spacetime correlations in the advecting velocity field
\cite{bj,pre}. Interestingly, neither diffusion nor advection when
acting alone can produce non-Gaussian tails in the PDF's. In our case,
the non-Gaussian tails would appear in the PDF of vortex displacements
during a fixed time interval $\Delta t$ (vortex flights).  In order to
test this possibility without detailed knowledge of the velocity field
${\bf V}({\bf r},t)$, we study the thermal motion of an isolated vortex
numerically as described below.

\section{Numerical Analysis and Results}

In this section we describe our numerical methods used to solve
Eqs.\ (\ref{dsnls}) under the constraint that only a single vortex be
present in the simulation volume. An {\it unpaired} vortex is forced
into the sample by coupling the field $\psi({\bf r},t)$ to a static
external gauge potential ${\bf A}({\bf r})$, replacing $\nabla$ in Eq.\ (6) by
$\nabla + i{\bf A}({\bf r})$. If we choose ${\bf A}({\bf r}) = -(2\pi/L^2)y\hat{x}$, there
will be on average one vortex per area $L\times L$ in an infinite 2D
sample. The introduction of ${\bf A}({\bf r})$ is a mathematical artifice and
does not represent a dynamical field. In the language of
superconductivity, ${\bf A}({\bf r})$ represents an external uniform magnetic
field that completely penetrates the sample without screening, and
vanishes in the thermodynamic limit $L\to\infty$.
This is a realistic approximation to a superconducting thin film geometry,
as in this case the external magnetic flux must fully penetrate the
superconductor.
Screening is only manifested as an inhomogeneity of the flux density
threading the plane of the superconductor. The characteristic length scale
for these variations is\cite{dG}
$\Lambda=\lambda^2/d$, where $\lambda$ is the bulk magnetic
penetration length of the superconducting material,
and $d$ is the thickness of the film. We study the system in
the limit where $\Lambda$ is much larger than any other length of the
problem.  In this limit the flux density threading the superconductor
is uniform and screening can be neglected.
One might also
worry that the derivation leading up to Eq.\ (18), which was done
explicitly for the case of the neutral superfluid, no longer applies
when ${\bf A}({\bf r})$ is introduced. While it is a fact that the
massless Goldstone boson in a superluid is in general replaced by a {\it
massive} vector boson in a superconductor, the mass of this boson
vanishes in the thin film limit considered,
and the gapless phonon spectrum is restored, as
pointed out by \v{S}im\'{a}nek.\cite{simanek}
By the same token the
electromagnetic vortex mass
vanishes in our thin film limit, as the moving vortex does
not generate any electric field.

We choose our system to be a square, $S=\langle 0,L)\times\langle 0,L)$.  We
wish to employ a variant of periodic boundary conditions so that the
entire $xy$ plane consists of physically equivalent copies of the
system at all times; by this we mean that the local superfluid density
$|\psi({\bf r},t)|^2$, and the supercurrent density
\begin{equation}
{\bf J}({\bf r},t) =
\frac{i}{2}(\psi\nabla\psi^*-\psi^*\nabla\psi)+|\psi|^2{\bf A}
\label{cond}
\end{equation}
should both be periodic functions with periodicity $L$ in both spatial
directions.  To this end we impose quasi-periodic boundary conditions
on $\psi({\bf r},t)$ (${\bf r}\in S$, $n_y$ and $n_y$ arbitrary integers):
\begin{equation}
\psi({\bf r}+n_xL\hat{x}+n_yL\hat{y},t)=\psi({\bf r},t)
\exp[in_x\theta_x+in_y(2\pi x/L+\theta_y)].
\end{equation}
The angles $\theta_x$ and $\theta_y$ are as yet arbitrary and may be
functions of time.  Upon discretization of Eqs.\ (\ref{dsnls}) we
observe that any given time dependence of $\theta_x$ and $\theta_y$
will in general affect the motion of $\psi({\bf r},t)$.  Yet $\theta_x$ and
$\theta_y$ are not dynamical degrees of freedom in the sense that
their motion is not determined by Eqs.\ (\ref{dsnls}). Therefore, we
are left with a problem of choosing the ``right'' time dependence for
$\theta_x$ and $\theta_y$. Fortunately the physical interpretation of
these variables is straightforward. In order to find it, we rewrite
the stochastic field equation in (\ref{dsnls}) in the form
\begin{equation}
(i/\gamma-1) \frac{\partial}{\partial t} \psi=
\frac{\delta H[\psi,\psi^*,{\bf A}]}{\delta \psi^*} +\epsilon\eta({\bf r},t),
\end{equation}
where
\begin{equation}
H[\psi,\psi^*,{\bf A}]\equiv \int
d^2{\bf r}\,\left[\frac{1}{2}|(\nabla+i{\bf A})\psi|^2
-|\psi|^2+\frac{1}{2}|\psi|^4\right]
\end{equation}
is the Hamiltonian of the system (the condensate) disconnected from
the thermal reservoir.  One can show that
\begin{equation}
\frac{\partial H}{\partial \theta_x}=
\hat{x}\cdot\int_0^Ldy\,\left.\frac{\delta H}{\delta {\bf A}({\bf r})}\right|_{x=0}
\equiv \hat{x}\cdot\int_0^Ldy\,\left.{\bf J}_s({\bf r},t)\right|_{x=0}.
\end{equation}
The first part of this equation is an algebraic equality; in the
second part we used the definition of the supercurrent, whose explicit
form appears in Eq.\ (\ref{cond}). Therefore, the first derivative of
the Hamiltonian with respect to $\theta_x$ is equal to the net
supercurrent passing through the edge of the sample perpendicular to
the $x$ axis. A completely analogous statement holds for the other
spatial direction as well.

If one is only interested in the thermal motion of the vortex in the
absence of external forces---such as, for example, the Magnus force
created by macroscopic superflow---the desired constraints on angles
$\theta_x$ and $\theta_y$ are
\begin{equation}
\frac{\partial H}{\partial \theta_x}=\frac{\partial H}{\partial
\theta_y}=0. \label{zfe}
\end{equation}
These constraints define a {\it zero flow ensemble} \cite{dpa},
relevant for experiments with no externally imposed driving flow.  The
angles $\theta_x$ and $\theta_y$, though not dynamical variables, have
their values fixed by two dynamical constraints, (\ref{zfe}), at all
times.

As a result of the quasi-periodic boundary conditions, the isolated
vortex is part of an infinite (always perfectly square) vortex
lattice, with lattice constant $L$.  The physical cutoff $\Lambda_c$
needed to evaluate the vortex mass is estimated as $\Lambda_c \simeq
L/\sqrt{\pi}$ (this is when a square of side $L$ and a circle of
radius $\Lambda_c$ have equal areas).  Our numerical approach consists
of simultaneously solving the lattice version of Eqs.\ (\ref{dsnls})
on the spatial domain $S$, and the two dynamical constraints
(\ref{zfe}).  The numerical grid spacing $\delta$ must resolve both
the vortex core size $\xi$, which is the unit of length, and the
phonon correlation length $\xi_p$, which is the characteristic length
scale of the four-point correlator $\langle
|\psi({\bf r},t)|^2|\psi({\bf r}',t)|^2\rangle$, and depends on temperature.  We
have found it adequate to set $\delta =0.15$ throughout the
temperature range studied.  Numerical stability of the second-order
stochastic Runge-Kutta method \cite{rk2} used to integrate
Eqs.\ (\ref{dsnls}) then allows a time step of $\tau=0.005$.  We study
sample sizes $L = 6$, 12, and 24, and temperatures $1/\beta = 0.0125$,
0.025, 0.05, and 0.1.

The vortex is clearly identifiable and no thermal vortex-antivortex
pairs are generated during the numerical runs. These ranged in
duration from $6\times 10^4$ to $6\times 10^5$, in dimensionless units
(up to $1.2\times 10^8$ time steps). The characteristic equilibration
time was $\sim 500$.  The vortex position is taken to be the
coordinate of the lower left-hand corner of the plaquette where it is
detected with the discretized version of the integral (\ref{vortpos}).

Fig.\ 1 shows $\langle [\Delta {\bf R}(t)]^2 \rangle$ for three
different $L$ and the same $\beta$.  It is apparent from the data
offsets that vortex motion is not purely diffusive. On the other hand,
the calculated diffusion coefficient $\sigma$ correctly approximates
the slopes of the numerical data at long times. Therefore we conclude
that the vortex mass, measured through $\sigma$, displays a
logarithmic divergence with the system size $L$.  The issue of the
value of the vortex inertial mass has been the subject of much
controversy in the literature \cite{nat,dl}.  Our results show that
even in the presence of random flow the logarithmic divergence of the
classical vortex mass persists. This stands in agreement with an
earlier result \cite{arovas} in the presence of coherent flow. It
would be interesting to investigate whether this divergence persists
in the presence of multiple vortices or in the limit of low damping.

An independent measure of diffusion is obtained through the concept of
the mean first passage time $\bar{t}(\ell)$ over the boundaries of the
square box of size $2\ell\times 2\ell$.  For pure diffusion one has
$\bar{t}(\ell) = 0.295\,\ell^2/\sigma$ \cite{sl}.  Fig.\ 2 reiterates
that the vortex motion is diffusive at long times. Deviations from
pure diffusion are manifest at short times that correspond to vortex
flights of ${\cal O}(1)$ and less.  The evidence from Fig.\ 2 also shows
that $\sigma$ depends linearly on temperature, as expected from
Eq.\ (\ref{diffcon}). This is a useful check on the quality of our
numerics---the vortex is not subject to spurious pinning due to the
discreteness of the underlying lattice. The flattening-out of the data
for large $\ell$ is due to finite sampling.

Finally in Fig.\ 3 we present the PDF of vortex flights. This
distribution has exponential tails (straight lines on a logarithmic
scale as in Fig. 3), in addition to the expected Gaussian profile
typical of Brownian motion. Our numerical results show that these
tails persist for all system sizes and $\beta$ studied, and the slope
decreases as the vortex mass $M(\Lambda_c)$ is increased (which
corresponds to a volume increase). The presence of exponential tails
has been linked to nontrivial time correlations in the advecting
velocity field \cite{bj}. The presence of finite correlation times can
be either a feature of a stochastic velocity field or can arise from
deterministic, chaotic motion (e.g., intermittency). In our case, we
expect the exponential tails to arise from the time correlations of
overdamped phonon modes which are allowed to exist in the fluctuation
spectrum of our theory. This aspect is presently under investigation.

\section{Discussion and Conclusions}

The particular models of vortex motion of the type described in the
Introduction typically ignore the presence of hydrodynamic modes.  In
contrast, real physical systems, in particular those described by
stochastic field theories, necessarily include these
effects. Consequently, as demonstrated here both analytically and
numerically, the motion of a single vortex must be viewed as diffusive
and randomly advected at the same time. The Fokker-Planck equation
(\ref{fpe}) makes explicit the formal analogy between single vortex
fluctuations and the dynamics of a passive scalar advected to random
fluid flow. The quantity $c({\bf r},t)$ is the ensemble-averaged
probability density of finding a vortex at ${\bf r}$ at an instant $t$.
The ensemble averaging is done over the fast degrees of freedom of the
heat bath. The slow, hydrodynamic, degrees of freedom of the
superfluid---the phonons---remain implicitly present in the velocity
field ${\bf V}({\bf r},t)$. This equation presents a picture of the vortex
as a diffusing massive particle subject to random advection. The
advection is apparently due to phonons, which propagate through the
vortex core randomly. As a result, the probability of large-scale
vortex flights is increased.

In turbulent fluid flow experiments with a dye as the passive scalar
it is possible to measure the dye concentration, analogous to our
$c({\bf r},t)$, directly, because individual dye particles do not interact
among themselves (apart from contact interaction).  Vortices, however,
interact logarithmically, so one cannot confirm the exponential tails
in $c({\bf r},t)$ directly by studying motion of a large aggregate of
vortices in a sample. Instead, one must study the motion of a single
vortex, and measure the statistics of vortex flights.

Examples where thermal vortex motion plays a key role include theories
of thermal depinning, flux creep, and dilute vortex lattice melting in
high-$T_c$ superconductors. Our current understanding of these
phenomena relies exclusively on the Brownian particle description.  It
would be interesting to re-examine these problems in the light of the
advected scalar aspect of vortex motion. An exciting new possibility
is a direct experimental observation of vortex motion in
two-dimensional high-$T_c$ superconductors by fast scanning tunneling
microscopy \cite{tak}. With this method the overwhelmingly more
frequent occurrence of large vortex flights in comparison with that of
a simple Brownian particle could be put to direct test.

\section{Acknowledgments}
We thank A.R. Bishop, D.P. Arovas, T. Hwa, and G. Lythe for useful
discussions.  Simulations were carried out at the Advanced Computing
Laboratory (ACL), Los Alamos National Laboratory, the Department of
Physics, Ohio State University, and at the National Energy Research
Scientific Computing Center (NERSC), Lawrence Berkeley National
Laboratory.