\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.


