text:
# Solutions of Gross–Pitaevskii equations beyond the hydrodynamic approximation: Application to the vortex problem.

## I Introduction

Since the remarkable experimental realization  of the Bose-Einstein condensation (BEC) there has been an explosion of the experimental and theoretical activity devoted to the physics of dilute ultracold bosonic gases .

From the theoretical point of view, many of current condensate systems are cold enough that the mean field theories can be applied to describe the ground state of the system . The most classical model arising in this context for the order parameter of the system is the the Gross-Pitaevskii (GP) equation. In fact, many relevant features of the condensate such as free expansion , low energy excitations , collapse , dynamics of multicomponent systems , vortex generation and dynamics , and many others, can be understood within the framework of this “simplest” mathematical model.

Mathematically, the GP equation, to be written explicitly below, is an equation of Nonlinear Schrödinger type. These equations have been studied extensively both in the physical and mathematical literatures  since they provide an universal model for the study of the dynamics of envelope waves. Some fields where these equations arise are plasma physics , fiber and integrated optics , wave propagation in Kerr media , water waves , biomolecule dynamics  and fundamentals of quantum mechanics . One of the distinctive features of the equation as it appears in BEC problems is the presence of an external potential -the trap-, which affects essentialy the spectrum of the “linear part”, i.e. of the underlying Schrödinger equation.

The mathematical analysis of the GP equation in the situations common in BEC is very complicated and this is why most theoretical studies deal with either numerical simulations or some kind of approximations such as the Thomas–Fermi, the hydrodynamic approach or simply reduce the number of degrees of freedom by using the variational approach based on trial functions. To be precise let us write explicitly the equation for the order parameter, whose form is

$$i\mathrm{}\frac{\mathrm{\Psi }}{t}=\frac{\mathrm{}^2}{2m}^2\mathrm{\Psi }+V_{ext}(𝐫)\mathrm{\Psi }+g|\mathrm{\Psi }|^2\mathrm{\Psi },$$
(1)
where the external potential can be well approximated by

$$V_{ext}=\frac{1}{2}m\left(\omega _x^2x^2+\omega _y^2y^2+\omega _z^2z^2\right),$$
(2)
and $`g=4\pi \mathrm{}^2a_s/m`$.

The simplest approximation that can be done to obtain information on the ground state is the so called Thomas–Fermi approximation, which corresponds to neglecting the kinetic energy term under the assumption that it is much less than the self-interaction energy. In the cases when the trap is smooth enough this approximation is not bad although it is not consistent, since the assumption of zero kinetic energy provides a profile that has an infinite value of that quantity.

The following level of approximation, the hydrodynamic approach, can be explained by introducing new variables $`n(𝐫,t)`$ and $`𝐯(𝐫,t)`$ through the relations

$`\mathrm{\Psi }(𝐫,t)`$ $`=`$ $`\sqrt{n(𝐫,t)}e^{iS(𝐫,t)},`$ (4)
$`𝐯(𝐫,t)`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{m}}S(𝐫,t).`$ (5)
Then, one arrives to the system of equations

$$_tn+(𝐯n)=0,$$
(6)
$$m_t𝐯+\left(V_{ext}+gn\frac{\mathrm{}^2}{2m\sqrt{n}}^2\sqrt{n}+\frac{mv^2}{2}\right)=0.$$
(7)
The conventional hydrodynamic approximation corresponds to neglecting the term proportional to $`^2\sqrt{n}`$ under the assumption that the density is a slowly varying function . For a spherically symmetric potential ($`\omega _x=\omega _y=\omega _z=\omega `$) this can be done in the limit $`Na_s/a_01`$, where $`a_s>0`$ is the scattering length and $`a_0=\sqrt{\mathrm{}/m\omega }`$ is the harmonic oscillator length, i.e. the trap size.

To the best of authors knowledge the mathematical limits of validity of the above mentioned approximations have not been established (which would imply study of the next orders of the expansion or making appropriate rigorous estimates), so far. In particular, using the estimate $`n10^{12}`$ cm<sup>-3</sup> one finds that for $`a_s10`$ nm the term with $`^2\sqrt{n}/\sqrt{n}`$ is much less than $`gn`$ for excitations whose characteristic length $`\mathrm{\Lambda }`$ is of order of $`10^3`$ cm. Such a $`\mathrm{\Lambda }`$ is comparable with the characteristic size $`a_0`$ of the trapping potential. Thus the question about analytical description of the excitations which vary on scales much less than $`a_0`$ appears. It is important to notice that in those cases the terms $`\frac{\mathrm{}^2}{m}\frac{^2\sqrt{n}}{\sqrt{n}}`$ and $`gn`$ become comparable with each other and have to be treated in equivalent manner. Using the wave terminology this situation can be characterized as balance of dispersive and nonlinear effects, which should lead to stable spatially localized excitations.

It is our intention in this paper to describe an analytical tool, the multiscale approximation, which is able to provide approximate solutions of the problems under consideration in explicit analytical form. In less rigorous terms this technique has been applied in the context of BEC in Refs.  to problems with negative scattering length. The application to positive scattering length presents some new characteristics as we will see later (in particular it implies different scaling order). Specifically, we will consider a BEC in a parabolic cilindrically or spherically symmetric confining potential and obtain a set of background solutions in the form of black holes or vortices whose characteristic spatial scale is much larger than the range of inter-atomic interactions but smaller than the trap size. Working in this way we will be able not only to obtain an approximate analytical shape for the vortex solutions but also to provide results on their dynamical stability, a point which is still of great interest after the many studies devoted to it .

The outline of the paper is as follows. First, in Sec. II we pose the problem in precise terms. As a second step, we apply the technique of multiscale expansion to our particular problem in Sec. III. In Sec. IV we obtain a particular set of solutions and discuss their stability. Finally, in Sec. V, we summarize our conclusions.

## II Statement of the problem

We start with the three-dimensional GP equation in the form

$$i\mathrm{}\frac{\mathrm{\Psi }}{t}=\frac{\mathrm{}^2}{2m}^2\mathrm{\Psi }+V_{ext}(𝐫)\mathrm{\Psi }+J_\lambda (𝐫)\mathrm{\Psi }.$$
(8)
We will preserve the integral form of the interaction term \[c.f. Eq. (1)\]

$`J_\lambda (𝐫)=g{\displaystyle V_\lambda (𝐫,𝐫^{})|\mathrm{\Psi }(𝐫^{},t)|^2𝑑𝐫^{}},`$ (9)
since there are situations where its role could be important  and it does not complicate the formalism. Here $`V_\lambda (𝐫,𝐫^{})`$ describes two-body interactions and $`\lambda `$ is the effective interaction length (mathematically, the size of the kernel).

The wave function is normalized to the total number of particles $`N`$

$$|\mathrm{\Psi }(𝐫)|^2𝑑𝐫=N.$$
(10)
It is reasonable to assume that the effective interaction length $`\lambda `$ parametrizes the kernel so that

$$\underset{\lambda 0}{lim}V_\lambda (𝐫,𝐫^{})=\delta (𝐫𝐫^{}).$$
(11)
The last limit corresponds to the usual case of local interactions and gives rise to the $`|\mathrm{\Psi }|^2`$ nonlinear term in the GP equation \[i.e. (8) transforms into (1)\]. Bearing in mind this last fact we can narrow the class of the kernels under consideration by requiring

$$V_\lambda (𝐫,𝐫^{})\frac{1}{\lambda ^3}U\left(\frac{|𝐫𝐫^{}|}{\lambda }\right),$$
(12)
where

$$U(|\stackrel{~}{𝐫}\stackrel{~}{𝐫}^{}|)𝑑\stackrel{~}{𝐫}^{}=1.$$
(13)
and $`\stackrel{~}{𝐫}=𝐫/\lambda `$. As it is clear from the definition of $`\lambda `$, it must be of order of the scattering length:

$$\frac{\lambda }{a_s}=O(1)$$
(14)
In what follows, however we keep $`\lambda `$ as a parameter.

We will concentrate on a BEC in a trap with either cilindrical (with $`z`$-axis being the symmetry axis) or spherical symmetry, i. e. $`\omega _\alpha =\omega `$, where $`\alpha =x,y`$ in the former case and $`\alpha =x,y,z`$ in the last case. We define the characteristic scales of the condensate in the trasversal, $`a_0=\sqrt{\mathrm{}/(m\omega )}`$, and longitudinal, $`a_z=\sqrt{\mathrm{}/(m\omega _z)}`$, directions. In the case of spherical symmetry these scales coincide, $`a_z=a_0`$.

In order to provide a formal description of the system it is convenient to scale out the GP equation. To this end we introduce the dimensionless time and wavefuntion through

$`\tau `$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{2m\lambda ^2}}t,`$ (16)
$`\psi (\stackrel{~}{𝐫},\tau )`$ $`=`$ $`2\lambda \sqrt{{\displaystyle \frac{2\pi a_s}{N}}}\mathrm{\Psi }(𝐫,t).`$ (17)
Then the renormalized GP equation takes the form

$$i\frac{\psi }{\tau }=\mathrm{\Delta }\psi +\stackrel{~}{V}_{ext}(\stackrel{~}{𝐫})\psi +\stackrel{~}{J}(\stackrel{~}{𝐫})\psi $$
(18)
where

$`\stackrel{~}{J}(\stackrel{~}{𝐫})={\displaystyle U(\stackrel{~}{𝐫}\stackrel{~}{𝐫}^{})|\psi (\stackrel{~}{𝐫}^{},\tau )|^2𝑑\stackrel{~}{𝐫}^{}},`$ (19)
$`\stackrel{~}{V}_{ext}(\stackrel{~}{𝐫})=\{\begin{array}{cc}(\mathrm{\Omega }^2\stackrel{~}{𝐫}_{}^2+\mathrm{\Omega }_z^2\stackrel{~}{z}^2)\hfill & \text{cilindrical symmetry}\hfill \\ \mathrm{\Omega }^2\stackrel{~}{𝐫}^2\hfill & \text{spherical symmetry}\hfill \end{array}`$ (22)
$`\stackrel{~}{𝐫}_{}=(\stackrel{~}{x},\stackrel{~}{y})`$, $`\mathrm{\Omega }=\lambda ^2/a_0^2`$ and $`\mathrm{\Omega }_z=\lambda ^2/a_z^2`$.

Having scaled out the kernel which describes the nonlocal interactions we have to impose the condition that the characteristic scale of the kernel variation is of order one, i. e.

$`\left|{\displaystyle \frac{1}{U(\stackrel{~}{r})}}{\displaystyle \frac{U(\stackrel{~}{r})}{\stackrel{~}{r}}}\right|=O(1).`$

where $`\stackrel{~}{r}=|\stackrel{~}{𝐫}|`$ and $`\stackrel{~}{𝐫}`$ must be substituted by $`\stackrel{~}{𝐫}_{}`$ in the case of the cilindrical symmetry.

## III Multiscale analysis.

The multiscale technique is a mathematical technique for the analysis of problems where there are different spatial scales involved in the solution. The problem at hand possesses three important scales: the scattering length $`a_s`$, the trap size $`a_0`$ and the healing length $`\xi =(8\pi na_s)^{1/2}`$. We will concentrate on the analysis of excitations of the superfluid corresponding to the scale of the healing length, which occur on top of a background with spatial scale of the order of the trap size. In other words our theory will be valid for cases when the above three parameters are related as follows

$$a_s\xi a_0.$$
(23)
For the first step we recall the properties of the linear spectral problem

$$L_d|n,l,m=_{n,l,m}|n,l,m,$$
(24)
where $`L_d=_d+_d`$. The operator $`_d`$ describes the radial dependence

$$_d=\frac{1}{\rho _0^{d1}}\frac{}{\rho }\rho _0^{d1}\frac{}{\rho _0}+\mathrm{\Omega }^2\rho _0^2,$$
(25)
and

$`_2`$ $`=`$ $`{\displaystyle \frac{^2}{z^2}}+{\displaystyle \frac{1}{\rho _0^2}}{\displaystyle \frac{^2}{\varphi ^2}},`$ (27)
$`_3`$ $`=`$ $`{\displaystyle \frac{1}{\rho _0^2\mathrm{sin}\theta }}{\displaystyle \frac{}{\theta }}\mathrm{sin}\theta {\displaystyle \frac{}{\theta }}+{\displaystyle \frac{1}{\rho _0^2}}{\displaystyle \frac{^2}{\varphi ^2}}.`$ (28)
The indexes $`n`$ and $`m`$ refer to radial and magnetic quantum numbers while $`l`$ is the azimutal (in the spherical case) or longitudinal (in the cilindrical case) quantum number (see Appendix). Hereafter $`d=2`$ and $`d=3`$ are associated with cilindrical and spherical symmetries.

For the next consideration it is convenient to single out the radial part of the eigenfunctions and represent

$`|n,l,m=\{\begin{array}{cc}\xi _{n,m}(\rho _0)\zeta _{l,m}(z,\varphi )\hfill & \text{cilindrical case}\hfill \\ \xi _{n,l}(\rho _0)\zeta _{l,m}(\theta ,\varphi )\hfill & \text{spherical case}\hfill \end{array}`$ (31)
To shorten notations in what follows we drop the second subindex of $`\xi `$ in the cases at hand: $`\xi _{n,}\xi _n`$, keeping in mind that $`n`$ is a radial (rather than total) quantum number.

The operator $`_d`$ is considered in the space of functions satisfying the conditions

$`|\xi _n(0)|<\mathrm{},\underset{\rho _0\mathrm{}}{lim}\rho _0\xi _n(\rho _0)=0.`$

The inner product in this space is defined by

$`(\xi _n,\xi _n^{})={\displaystyle _0^{\mathrm{}}}\overline{\xi }_n(\rho )\xi _n^{}(\rho )\rho ^{d1}𝑑\rho .`$ (32)
Then the operator $``$ is Hermitiam.

The problem at hand is directly related to the well known eigenvalue problem for the linear oscillator in two and three dimensional parabolic potential . Its spectrum is discrete

$$_{n,l,m}=\{\begin{array}{cc}4\mathrm{\Omega }\left(n+\frac{1}{2}+\frac{m}{2}\right)+2\mathrm{\Omega }_z\left(l+\frac{1}{2}\right)\hfill & \text{cilindical case}\hfill \\ 4\mathrm{\Omega }\left(n+\frac{3}{4}+\frac{l}{2}\right)\hfill & \text{spherical case}\hfill \end{array}$$
(33)
(both cases naturally coincide after substitution $`\mathrm{\Omega }_z\mathrm{\Omega }`$ and $`l+ml`$ in the first formula). The orthormalized eigenfunctions of $`_d`$ are given in the Appendix.

In accordance with the multiscale method, which is a theory for weakly nonlinear problems, the solution is searched in the form of a weakly modulated linear mode. This is a main feature of the multiscale method, which means that the theory is applicable to small condensates. However, in many problems the results coming out from the multiscale expansions remain valid even for the large nonlinearity limit. Respectively we introduce a formal small parameter $`ϵ1`$ and look for the solution of (18) in the form

$$\psi =ϵ\psi _1+ϵ^2\psi _2+\mathrm{}=\left(ϵ\xi _n^{(1)}+ϵ^2\xi _n^{(2)}+\mathrm{}\right)\zeta _{l,m}$$
(34)
Some peculiarities of the expansion introduced are to be discussed here. First, we use the expansion of only radial part of the solution. This means that the consideration will be restricted to symmetric solutions (the generalization is straightforward but cumbersome). Second, in what follows we will concentrate on modulation of the ground state which turns out to be flat compared with the excitation itself. This however is not in contradiction with the assumption about smooth (i.e. dependent on $`\rho _1`$) modulation of the underlying linear mode, depending on the rapid variable $`\rho _0`$, by weak nonlinearity. This is justified by the following reasons: (i) The background solution must possess the “right” assymptotic behavior at $`\rho \mathrm{}`$, which can be obtained only by taking into account the term $`\mathrm{\Omega }^2\rho ^2`$ at the “linear level”; (ii) The expansion must take into account the contribution of all linear eigenmodes, rather than only the leading one \[see (41) and discussion below\]; (iii) The linear problem with the parabolic trap is solved exactly.

The functions $`\xi _n^{(j)}`$ depend on the set of spatial coordinates $`\rho _j=ϵ^j\rho `$ and times $`\tau _j=ϵ^j\tau `$ ($`j=0,1,2,\mathrm{}`$): $`\xi _n^{(j)}\xi _n^{(j)}(\rho _0,\rho _1,\rho _2,\mathrm{};\tau _0,\tau _1,\tau _2,\mathrm{})`$. Hereafter $`\rho =|\stackrel{~}{𝐫}_{}|`$ in the cilindrical case and $`\rho =|\stackrel{~}{𝐫}|`$ in the spherical case. The variables $`\rho _j`$ and $`\tau _j`$ are regarded as independent: thus

$`{\displaystyle \frac{}{\tau }}`$ $`=`$ $`{\displaystyle \frac{}{\tau _0}}+ϵ{\displaystyle \frac{}{\tau _1}}+ϵ^2{\displaystyle \frac{}{\tau _2}}+O(ϵ^3)`$ (36)
$`{\displaystyle \frac{1}{\rho ^{d1}}}{\displaystyle \frac{}{\rho }}\rho ^{d1}{\displaystyle \frac{}{\rho }}`$ $`=`$ $`{\displaystyle \frac{1}{\rho _0^{d1}}}{\displaystyle \frac{}{\rho _0}}\rho _0^{d1}{\displaystyle \frac{}{\rho _0}}+ϵ\left(2{\displaystyle \frac{^2}{\rho _0\rho _1}}+{\displaystyle \frac{d1}{\rho _0}}{\displaystyle \frac{}{\rho _1}}\right)+O(ϵ^2)`$ (37)
Before going into details of the multiscale expansion we have to clarify the physical meaning of the parameter $`ϵ`$. We are interested in excitations against a background. Then in terms of the dimensionless wave function $`\psi _1`$ the normalization condition (10) leads to the estimate $`N\frac{ϵ^2V}{8\pi a_s\lambda ^2}`$ where $`V=a_0^da_z^{3d}`$ is the volume of the condensate which is limited by the confining potential. These estimates imply the normalization condition

$$|\psi _1(\stackrel{~}{𝐫},\tau )|^2𝑑\stackrel{~}{𝐫}=\frac{a_0^da_z^{3d}}{\lambda ^d}$$
(38)
(notice that in the case of cilindrical symmetry, taking $`\lambda =a_s`$ one obtains that the new wave function is normalized to one). In other words the effective small parameter of the multiscale expansion can be identified as

$$ϵ=O\left(\sqrt{na_s^3}\lambda /a_s\right)=O\left(\sqrt{na_s^3}\right)$$
(39)
where $`n=N/V`$ is the density of the condensate and (14) is used. Taking into account that typically $`na_s^3=O(10^4÷10^6)`$ one concludes that $`ϵ`$ ($`10^2÷10^3`$) can indeed be used for the multiscale expansion. We notice that the approximation of the local kernel, $`\lambda 0`$ is evidently a particular limit of the theory. Formally, however, it is more convenient to keep $`ϵ`$ as a small parameter without specification and substitute it by one in the final formulas.

We are interested in the evolution of an initially unperturbed ”linear” state \[say, $`(n_0,l_0,m_0)`$th one\]. Gathering the terms of the first order with respect to $`ϵ`$ expansion we conclude that $`\psi _1`$ can be represented in the form ($`_0=_{n_0,l_0,m_0}`$)

$$\xi _0^{(1)}=A(\rho _1,\mathrm{};\tau _1,\mathrm{})e^{i_0\tau _0}|n_0,l_0,m_0.$$
(40)
Hereafter the notation $`(\rho _1,\mathrm{};\tau _1,\mathrm{})`$ is used in order to indicate the dependence on all spatial and temporal variables, which are slower than those written explicitely \[i.e. for example in (40) the amplitude depends on $`\tau _1,\tau _2,\mathrm{}.`$\].

The solution of the equation appearing in the second order of $`ϵ`$ is represented in the form of the expansion

$$\psi _2=\underset{n,l,m}{\overset{}{}}B_n^{}(\rho _1,\mathrm{};\tau _1,\mathrm{})e^{i_0\tau _0}|n,l,m.$$
(41)
Hereafter a prime means that the sum is computed over all arguments such that $`nn_0`$, $`ll_0`$, and $`mm_0`$. The respective equation of the second order takes the form

$$\underset{n,l,m}{\overset{}{}}(_0_{n,l,m})B_{n,l,m}|n,l,m+\left[i\frac{A_n}{\tau _1}+\left(2\frac{}{\rho _0}+\frac{d1}{\rho _0}\right)\frac{A_n}{\rho _1}\right]|n_0,l_0,m_0=0$$
(42)
Applying $`n_0,l_0,m_0|`$ to this equation (see Appendix) one arrives at the relation $`A/\tau _1=0`$. Here we use the property

$`(\xi _n,\widehat{p}\xi _n)=i{\displaystyle _0^{\mathrm{}}}\xi _n(\rho _0)\left(2{\displaystyle \frac{}{\rho _0}}+{\displaystyle \frac{d1}{\rho _0}}\right)\xi _n(\rho _0)\rho _0^{d1}𝑑\rho _0=0`$

with

$$\widehat{p}=i\left(2\frac{}{\rho _0}+\frac{d1}{\rho _0}\right)$$
(43)
being the operator of the radial component of the (dimensionless) linear momentum. Then one concludes that $`A`$ does not depend on $`\tau _1`$, $`AA(\rho _1,\mathrm{};\tau _2,\mathrm{})`$.

The coefficient $`B_{n_0,l_0,m_0}`$ is equal to zero identically \[otherwise it could be made zero by simple renormalization of initial conditions for (18)\]. In order to find other coefficients $`B_{n,l,m}`$ we compute the inner product with $`n_0,l_0,m_0|`$. Introducing matrix elements

$$p_{n^{}n}=(\xi _n,\widehat{p}\xi _n^{})$$
(44)
corresponding to the transition between the states $`n^{}`$ and $`n`$, originated by the radial component of the linear momentum considered as a perturbation, we obtain the coefficients of the second order term

$$B_{n,l,m}(\rho _1,\mathrm{};\tau _2,\mathrm{})=i\frac{p_{n_0n}}{_0_{n,l,m}}\frac{A}{\rho _1}$$
(45)
Finally, collecting all the terms of the third order of $`ϵ`$ and using the explicit form of the second order addendum found above we arrive at the dynamical equation for the amplitude $`A`$

$$i\frac{A}{\tau _2}+𝒟_0\frac{^2A}{\rho _1^2}\chi _{n_0n_0}|A|^2A=0$$
(46)
Here

$$𝒟_0=𝒟_{n_0,l_0,m_0}=1\underset{n,l,m}{\overset{}{}}\frac{|p_{n_0n}|^2}{_{n,l,m}_0}$$
(47)
is the effective dispersion and

$$\chi _{nn^{}}=\xi _n^2(\rho )\xi _n^{}^2(\rho ^{})U(|\stackrel{~}{𝐫}\stackrel{~}{𝐫}^{}|)𝑑\stackrel{~}{𝐫}𝑑\stackrel{~}{𝐫}^{}$$
(48)
is the effective nonlinearity. As a matter of fact the effective nonlinearity can be computed under the assumption that the characteristic scale of the ground state is much larger than the width of the kernel $`U(|\stackrel{~}{𝐫}|)`$ of the interaction potential. Taking into account (14), (23), and (A1), (A3) one concludes that this supposition is valid with great accuracy for several lowest eigenstates $`\xi _n(\rho )`$. Then $`U(|𝐫|)`$ in the integrand of (48) can be substituted by the delta function and thus $`\chi _{nn}=(\xi _n^2,\xi _n^2)`$.

In the explicit form of $`𝒟_0`$ one can recognize the second order addendum to the energy of an effective particle induced by the perturbation operator $`\widehat{p}`$. Thus the effective dispersion has a contribution of the probability of the energy transfer among the levels in the confining potential. Notice that the mentioned transition occur among levels with the same orbital and azimutal quantum numbers (in the case at hand they are equal to zero).

Equation (46) is considered on the semiline ($`\rho _10`$) subject to the boundary condition $`lim_{\rho _1\mathrm{}}A(\rho _1,\mathrm{};\tau _2,\mathrm{})=𝒜e^{i\stackrel{~}{\omega }\tau _2}`$ where the constant $`𝒜`$ can be interpreted as an amplitude of the fundamental state of the condensate: after all it must be found from the normalization condition. \[Notice that other boundary conditions, of less practical interest, are also compatible with the problem at hand\].

Let us now return to the discussion of the main scales of the problem. In the dimensionless variables the spatial scale of the amplitude $`A(\rho _1,\mathrm{};\tau _2,\mathrm{})`$ can be estimated as $`1/ϵ`$. On the other hand the problem possess another important dimensionless scale $`1/\sqrt{\mathrm{\Omega }}`$ which determines the wave function localization due to the confining potential. An interesting and physically relevant situation appears when $`ϵ\sqrt{\mathrm{\Omega }}`$, what in physical variables means $`8\pi Na_s/a_01`$ , a condition which is satisfied in current experimental setups. The last requirement can be understood equivalently as smallness of the healing length $`\xi `$ compared with the linear oscillator length $`\xi a_0`$. In this situation the excitation width is much smaller than the volume “available” for the condensate and in the case of positive scattering length one can obtain localized excitations against the background. This is the regime of validity of our theory as has been put in more mathematical terms in the first part of this section.

## IV Application: Vortex shape and stability

One of the simplest excitations of the ground state condensate one can construct have the form of static ”dark” soliton solutions of Eq. (46), i. e. “black holes” or vortices (when a local phase is added). Joining the different contributions from the multiscale method we obtain

$`\mathrm{\Psi }={\displaystyle \frac{1}{\sqrt{\pi }}}{\displaystyle \frac{a_s}{\xi }}\mathrm{exp}\left[i\left({\displaystyle \frac{2}{a_0^2}}+{\displaystyle \frac{1}{a_z^2}}+8\pi \sqrt{{\displaystyle \frac{n}{a_s}}}\right){\displaystyle \frac{\mathrm{}t}{2m}}\right]\mathrm{exp}\left({\displaystyle \frac{r_{}^2}{2a_0^2}}{\displaystyle \frac{z^2}{2a_z^2}}\right)\mathrm{tanh}\left({\displaystyle \frac{r_{}}{\sqrt{2𝒟_0}\xi }}\right)e^{im\varphi }+O(na_s^3).`$ (49)
This solution approximates the shape of a vortex line with axis of rotation $`z`$ in trap of transverse size $`a_0`$.

This is not the only solution one can construct, as an example we present also a solution corresponding to a “black hole” in three dimensions

$`\mathrm{\Psi }={\displaystyle \frac{1}{\pi ^{3/4}}}{\displaystyle \frac{a_s}{\xi }}\mathrm{exp}\left[i\left({\displaystyle \frac{3}{a_0^2}}+2^{7/2}\sqrt{{\displaystyle \frac{\pi n}{a_s}}}\right){\displaystyle \frac{\mathrm{}t}{2m}}\right]\mathrm{exp}\left({\displaystyle \frac{r^2}{2a_0^2}}\right)\mathrm{tanh}\left(\left({\displaystyle \frac{2}{\pi }}\right)^{1/4}{\displaystyle \frac{r}{\sqrt{2𝒟_0}\xi }}\right)+O(na_s^3)`$ (50)
where $`r=(x,y,z)`$.

In the above formulas we have neglected the difference between $`\lambda `$ and $`a_s`$ \[see (14)\]. They however still contain the quantity $`𝒟_0`$ to be estimated. To this end we employ (A1) and (A3) and compute (47) . The result is $`𝒟_00.586`$ and $`𝒟_00.728`$ for (49) and (50) respectively. Thus in the physical units the effective width of the black hole can be estimated to be of order of the healing lengths which is the expected result for the vortex size.

It must be said that the result naturally does not depend on $`\mathrm{\Omega }`$, a fact which is explained by taking into account that in the parabolic trap the energy levels are equidistant. Another point to be mentioned is that as a matter of fact we have computed the terms of order of $`O(na_s^3)`$. They are given by (41), (45) and are not represented explicitly because of cumbersome form. In other words the theory provides one with the BEC wavefunction computed with accuracy $`O((na_s^3)^{3/2})`$ to problems whose number of particles is not too large.

One of the advantages of the multiscale method is that it guarantees the dynamical stability of the solution up to a certain time scale. In order to explain this last point we notice that the method is based on eliminating the secular terms in the third order with respect to the small parameter. This in particular means that the characteristic time of evolution of instabilities of the solutions (if any) is given by $`ϵ^3\tau `$. The main implication is that the instabilities against the background solution (which is the vortex in the case at hand) cannot develop up to $`\tau ϵ^3`$. Recalling now (39) and returning to the physical variables we obtain that for the typical values of $`ϵ`$ $``$ 10<sup>-2</sup> any instability could develop only for times larger than 10 seconds, which is of the order of the life-time of the condensate. In fact, this prediction is consistent with the previous numerical observation  that the vortex is a stable object.

This result is relevant since a mathematical proof of the expected dynamical stability of a vortex in a trap (i.e. when dissipation is zero) is not available yet. However, the numerical simulations (which are also valid only for limited times) point that a vortex should be a dynamically stable object. The analysis of the linearized operator around the vortex done in Ref.  was not completely conclusive since there were some zero eigenvalues which could lead to instabilities mediated by high order terms in the linearized evolution operator. In this sense the multiscale technique adds a new argument in favor of long time stability of the vortex in the framework of the mean field modelization of the phenomenon.

We remark that the previous discussion applies to dynamical stability, i. e. stability of the solution under perturbations on the initial data. It was discussed on the previous works on vortex stability that the system is energetically unstable , while the question of its Lyapunov stability remains open. Our analysis concerns only the subject of dynamical stability and confirm the numerical observations by many authors that vortices in a trap are very long lived and could even be completely stable .

## V Conclusion

In the present paper we have developed the multiscale approach for obtaining explicit (approximate) solutions to the GP equation for condensates with sufficiently small number of particles which are characterized by fast spatial scales of the order of the healing length, and thus cannot be described within the framework of the conventional hydrodynamical approach. One of the advantages of the method is that it guarantees the stability of the solution as discussed in Sec. IV which allows us to put lower bounds on the time a vortex survives under possible dynamical instabilities described within the framework of the GP equation.

Another outcome of the theory developed above is that the solutions of the healing length size cannot be excited only as a background state. Some accompanion modes, corresponding to higher levels of the confining potential are excited simultaneously. All these contributions are included in the constants $`𝒟_0`$ appearing in the reduced equations (46).

One more advantage of the multiscale method is that, by reducing the problem to solving the effective one dimensional nonlinear Schrödinger equation, one gets a powerful tool for obtaining other more complicated solutions and even interactions among them, a point which will be the subject of further study.

The present paper provides a new analytical tool which we hope will be of interest for theoretical progresses in the analysis of BEC problems.

###### Acknowledgements.

We are in debt to G. Alfimov for many fruitful discussions. V. V. K. is grateful to University of Castilla-La Mancha for warm hospitality. V.M.P-G. has been supported by DGCYT under grant PB96-0534. The cooperative work has been supported through the bilateral program DGCyT-HP1999-019/Ação N<sup>o</sup> E-89/00.

## A Eigenfunctions

For the sake of convenience here we present the orthonormalized eigenfunctions. In the case of potential with cilindrical symmetry they are

$$\xi _{n,m}(\rho )=\sqrt{2\mathrm{\Omega }\frac{n!}{(n+m)!}}\left(\sqrt{\mathrm{\Omega }}\rho \right)^me^{\mathrm{\Omega }\rho ^2/2}L_n^{(m)}(\mathrm{\Omega }\rho ^2),$$
(A1)
$$\zeta _{l,m}(\varphi ,z)=\frac{\mathrm{\Omega }_z^{1/4}}{\pi ^{3/4}\sqrt{2^{l+1}l!}}e^{im\varphi }e^{\mathrm{\Omega }_zz^2/2}H_l\left(\sqrt{\mathrm{\Omega }}_zz\right).$$
(A2)
Here $`H_n()`$ is the Hermite polinomial and $`L_n^{(m)}()`$ is the Laguerre polinomial. When spherical symmetry is imposed the eigenfunctions are given by

$`\xi _{n,l}(\rho )=2\left[{\displaystyle \frac{2^{n+l}}{(2(n+l)+1)!!}}\mathrm{\Gamma }(n+1)\right]^{1/2}\left({\displaystyle \frac{\mathrm{\Omega }^3}{\pi }}\right)^{1/4}e^{\mathrm{\Omega }\rho ^2/2}(\sqrt{\mathrm{\Omega }}\rho )^lL_n^{(l+1/2)}(\mathrm{\Omega }r^2),`$ (A3)
$`\zeta _{l,m}(\theta ,\varphi )Y_l^m(\theta ,\varphi )=(1)^{m+|m|}i^l\left[{\displaystyle \frac{2l+1}{4\pi }}{\displaystyle \frac{(l|m|)!}{(l+|m|)!}}\right]^{1/2}P_l^{|m|}(\mathrm{cos}\theta )e^{im\varphi }.`$ (A4)
The inner product in this space is defined by

$$n,l,m|n^{},l^{},m^{}=_{\mathrm{}}^{\mathrm{}}𝑑z_0^{2\pi }𝑑\varphi _0^{\mathrm{}}𝑑\rho \rho \overline{\xi }_{n,m}(\rho )\xi _{n^{},m^{}}(\rho )\overline{\zeta }_{l,m}(z,\varphi )\zeta _{l^{},m^{}}(z,\varphi )$$
(A5)
in the case of cilindrical symmetry and by

$$n,l,m|n^{},l^{},m^{}=_\pi ^\pi 𝑑\theta _0^{2\pi }𝑑\varphi _0^{\mathrm{}}𝑑\rho \rho \overline{\xi }_{n,m}(\rho )\xi _{n^{},m^{}}(\rho )\overline{\zeta }_{l,m}(\theta ,\varphi )\zeta _{l^{},m^{}}(\theta ,\varphi )$$
(A6)
in the case of spherical simmetry


format: text

text:
# Strangeness in Lattice QCD

## 1 Introduction

Quantum Chromodynamics (QCD) is widely accepted as the theory describing the strong interaction. It is formulated in terms of quarks and gluons and contains relatively few free parameters, namely the gauge coupling $`g`$ and the masses $`m_\mathrm{u},m_\mathrm{d},m_\mathrm{s},\mathrm{}`$ of the various quarks. These parameters must be fixed by using experimental input: a typical example is the determination of the sum of quark masses $`(m_\mathrm{u}+m_\mathrm{s})`$, using the experimental value of the kaon mass, $`m_{\mathrm{K}^+}`$:

$$m_{\mathrm{K}^+}(m_\mathrm{u}+m_\mathrm{s}).$$
(1)
Clearly, in order to connect the two sides of this relation one has to “solve QCD” for the case at hand.

The QCD coupling constant $`g`$, which describes the coupling strength, depends on the momentum transfer $`q`$ between quarks and gluons. The well-known property of “asymptotic freedom” implies that $`g`$ decreases for large $`q`$. By contrast, the coupling becomes large when momentum transfers of the order of typical hadronic scales are considered, say, $`q1\mathrm{GeV}/c`$. As a consequence, perturbation theory in $`g^2`$ becomes an inadequate tool to study QCD at low energies. Hence, a non-perturbative treatment is required in order to understand hadronic properties on a quantitative level. The formulation of QCD on a (Euclidean) space-time lattice  provides such a non-perturbative method to compute the relations between experimentally accessible quantities and the parameters of QCD through numerical simulations.

Here I restrict myself to QCD at zero temperature, focusing on the spectrum of light hadrons and the determination of the strange quark mass from first principles. A more general overview of the status of lattice calculations can be obtained from the proceedings of the annual Conference on Lattice Field Theory  and other review talks .

In the next section I shall outline the basic concepts of the lattice formulation. Section 3 contains results from recent benchmark calculations of the light hadron spectrum. In section 4 recent determinations of the strange quark mass are described. Finally, section 5 contains some concluding remarks.

## 2 Basic concepts

The lattice formulation of QCD replaces the familiar continuous Minkowski space-time by a discretized, Euclidean version with finite volume $`L^3T`$. Points in space-time are separated by a finite distance $`a`$, the lattice spacing. The inverse lattice spacing, $`a^1`$, acts as an UV cutoff, which regularizes the infinities that are typically encountered in Quantum Field Theory. The quark and antiquark fields $`\psi (x),\overline{\psi }(x)`$ are associated with the sites of the lattice, whereas the gauge field is represented by the so-called link variable $`U_\mu (x)`$, which connects neighbouring lattice sites, and is an element of the gauge group SU(3). After choosing a suitable, gauge-invariant discretization of the QCD action

$$S[U,\overline{\psi },\psi ]=S_\mathrm{G}[U]+S_\mathrm{F}[U,\overline{\psi },\psi ],$$
(2)
where $`S_\mathrm{G}`$ is the lattice action of pure Yang–Mills theory and $`S_\mathrm{F}`$ denotes the lattice fermion action, one may define the expectation value of an observable $`\mathrm{\Omega }`$ as

$$\mathrm{\Omega }=\frac{1}{Z}D[U]D[\overline{\psi }]D[\psi ]\mathrm{\Omega }\mathrm{exp}\left(S_\mathrm{G}S_\mathrm{F}\right).$$
(3)
Here the normalization of the functional integral $`Z`$ is determined by requiring $`11=1`$. After performing the integration over fermionic fields the expression for the expectation value becomes

$$\mathrm{\Omega }=\frac{1}{Z}\underset{x,\mu }{}dU_\mu (x)\underset{\mathrm{f}}{}det\left(D+m_\mathrm{f}\right)\mathrm{\Omega }\mathrm{exp}(S_\mathrm{G}),$$
(4)
where $`D`$ is the lattice Dirac operator and $`m_\mathrm{f}`$ is the mass of quark flavour $`f`$. Thus, the discretization procedure has given a meaning to the functional integration over gauge fields, which reduces to an ordinary multiple-dimensional integral over group elements. The various steps leading to eq. (4) can be treated in a gauge-invariant manner, and hence the lattice formulation represents a regularization procedure that preserves the gauge invariance of QCD. Furthermore, the definition of physical observables does not rely on perturbation theory. Therefore eq. (4) forms the basis for a non-perturbative, stochastic evaluation of expectation values of physical observables using numerical simulations.

Despite an enormous increase in computer power, there remain several major difficulties that make realistic simulations of QCD a hard task. Perhaps the biggest challenge is the inclusion of dynamical quark effects: the evaluation of the fermionic determinant in eq. (4) in numerical simulations is still very costly, even on today’s massively parallel computers. In early simulations of lattice QCD the determinant was therefore set equal to 1, a choice whose physical interpretation corresponds to neglecting quark loops in the evaluation of $`\mathrm{\Omega }`$. Although this represents a rather drastic assumption about the influence of quark-induced quantum effects, the quenched approximation works surprisingly well, as we shall see later. However, it is clearly desirable to develop simulation algorithms that allow for a more efficient evaluation of the fermionic determinant. With present algorithms, the cost of simulating “full” QCD is roughly a thousand times higher than that of the corresponding quenched simulation.

An indirect consequence of using the quenched approximation is the so-called scale ambiguity. That is, the calibration of the lattice spacing in physical units, $`a^1[\mathrm{MeV}]`$, depends on the quantity $`Q`$, which is used to set the scale

$$a^1[\mathrm{MeV}]=\frac{Q[\mathrm{MeV}]}{(aQ)},Q=f_\pi ,m_N,m_\rho ,\mathrm{}$$
(5)
This ambiguity arises because different quantities $`Q`$ are affected by quark loops in different ways. Clearly, dimensionful quantities such as particle masses are directly affected by the scale ambiguity.

Another problem one has to address are lattice artefacts (cutoff effects). Let $`\mathrm{\Omega }`$ denote the mass of a hadron in units of, say, the nucleon mass. Then the expectation values on the lattice and in the continuum differ by corrections of order $`a^p`$:

$$\mathrm{\Omega }^{\mathrm{lat}}=\mathrm{\Omega }^{\mathrm{cont}}+O(a^p),$$
(6)
where the value of integer $`p`$ in the correction term depends on the chosen discretization of the QCD action. Values of $`a`$ that can currently be simulated lie in the range $`a0.2`$$`0.05\mathrm{fm}`$. The size of the correction term can in some cases be as large as $`20`$%, depending on the quantity and the chosen discretization. An extrapolation to the continuum limit, $`a0`$, is then required to obtain the desired result. Surely this extrapolation is much better controlled if the discretization avoids small values of $`p`$.

Finally there are restrictions on the quark masses $`m_\mathrm{f}`$ that can be simulated. In general the following inequalities should be satisfied in any simulation:

$$a\xi L,$$
(7)
where $`L`$ is the spatial extent of the lattice volume. The quantity $`\xi `$ denotes the correlation length of a typical hadronic state and serves as a measure of the quark mass. The inequality on the right places restrictions on the light quark masses that can be simulated: if those are too light one may suffer from finite-size effects, since $`\xi `$ becomes large. Typical spatial extensions of $`L1.5`$$`3\mathrm{fm}`$ imply that the physical pion mass cannot be reached. The left inequality restricts the masses of heavy quarks. Since $`a^12`$$`4\mathrm{GeV}`$, it is clear that relativistic $`b`$-quarks cannot be simulated. One therefore relies on extrapolations in $`m_\mathrm{f}`$ to connect to the physical $`u,d`$ and $`b`$ quarks.

## 3 Hadron spectroscopy

In many ways the calculation of the mass spectrum of light hadrons is a benchmark of lattice QCD. A comparison of the experimentally observed spectrum with the results from an “ideal” lattice simulation of QCD, in which dynamical quark effects, lattice artefacts and quark mass dependences are all sufficiently controlled, would represent a stringent test of QCD as the underlying theory of the strong interaction. Alternatively, such a comparison enables us to assess the inherent systematics of current simulations, notably the effects of dynamical quarks. In particular one may investigate the quality of the quenched approximation.

### 3.1 Quenched light hadron spectrum

Recently the CP-PACS Collaboration presented a precision calculation of the light hadron spectrum in quenched QCD , which superseded an earlier study by GF11 . The findings of CP-PACS are summarized in the plot shown in Fig. 1. The main features of the spectrum are well reproduced by the quenched lattice data. Bearing in mind that only two input quantities have been used, namely the mass of the $`\rho `$-meson to set the scale and the mass of a strange meson (either $`m_\mathrm{K}`$ or $`m_\varphi `$) to fix the mass of the strange quark, the fact that the masses of so many hadrons are predicted quite accurately represents a major achievement in the understanding of QCD. Nevertheless, one finds small but significant deviations from the experimentally observed spectrum. For instance, the ratio of the nucleon and $`\rho `$ masses is calculated as

$$m_\mathrm{N}/m_\rho =1.143\pm 0.033,$$
(8)
which is 6.7% (2.5$`\sigma `$) below the experimental value of 1.218. Similarly, vector-pseudoscalar mass splittings such as $`m_\mathrm{K}^{}m_\mathrm{K}`$ are too small by $`10`$$`16`$% (4–6$`\sigma `$), depending on whether $`m_\mathrm{K}`$ or $`m_\varphi `$ is used to fix the strange quark mass. This implies that, for the first time, a significant deviation between the quenched QCD spectrum and nature is detected. Thus the conclusion is that quenched QCD describes the light hadron spectrum at the level of 10%. However, it also shows that the quenched approximation works surprisingly well, since the discrepancy is fairly mild. This has important consequences for lattice predictions of some phenomenologically interesting quantities, for which one has to rely on the quenched approximation for some time.

Although the CP-PACS results represent a real benchmark in terms of statistics, parameter values and lattice volumes, further corroboration of these findings is required. Recent calculations employing different discretizations have largely confirmed the findings of : the MILC Collaboration  has used staggered fermions  and finds a value for $`m_\mathrm{N}/m_\rho `$ in the continuum limit of

$$m_\mathrm{N}/m_\rho =1.254\pm 0.018(\mathrm{stat})\pm 0.028(\mathrm{syst}).$$
(9)
This is in broad agreement with experiment, but the difference to the CP-PACS result amounts to only 2$`\sigma `$. A recent calculation by UKQCD using $`O(a)`$ improved Wilson fermions  finds $`m_\mathrm{N}/m_\rho =1.26_{14}^{+8}`$. Whereas the error is too large to detect a significant deviation, UKQCD’s results for the spectrum also indicate that the quenched light hadron spectrum agrees with experiment at the level of 10%.

### 3.2 Beyond the quenched approximation

An obvious question is whether sea quark effects can account for the observed deviation of the quenched light hadron spectrum from experiment. Recently several collaborations have studied the light hadron spectrum using $`N_\mathrm{f}=2`$ flavours of dynamical quarks, which are identified with the physical up and down quarks . However, the masses of the dynamical quarks that can be simulated are quite large. A measure for the difference between the light quark masses found in nature and those used in dynamical simulations is provided by the pseudoscalar-to-vector mass ratio, $`m_{\mathrm{PS}}/m_\mathrm{V}`$. Current simulations are typically performed for $`m_{\mathrm{PS}}/m_\mathrm{V}=0.6`$–0.8, whereas the physical value is $`m_{\mathrm{PS}}/m_\mathrm{V}=m_\pi /m_\rho =0.169`$. Therefore one relies on extrapolations in the sea quark mass to make contact with the physical case. In addition, it has not been possible to simulate a third dynamical quark with $`m_\mathrm{f}m_\mathrm{s}`$ efficiently using known algorithms. Therefore, in most dynamical calculations an unphysical number of sea quarks has been used.

Despite these shortcomings one can perform systematic studies of the influence of dynamical quarks on the light hadron spectrum. CP-PACS  have studied the continuum limit of hadron masses computed for $`N_\mathrm{f}=2`$ and compared it with the results of the quenched light hadron spectrum discussed before. They find that the discrepancy with experiment is reduced from 6.5% to 1.4% for $`m_\varphi `$ computed with $`N_\mathrm{f}=2`$. For $`m_\mathrm{K}^{}`$ the gap decreases from 4.4% to 1.0%. This demonstrates clearly that the experimentally observed spectrum is reproduced more closely when sea quarks are “switched on”. Furthermore, results for the vector-pseudoscalar mass splitting reported by UKQCD  also show that lattice data for this quantity approach the experimental value as the sea quark mass is decreased.

Of course, the small remaining differences between QCD with $`N_\mathrm{f}=2`$ and experiment have to be explained. It is reasonable to assume that the unphysical value of $`N_\mathrm{f}`$ will have some influence. Further studies are also required to decide whether there is yet sufficient control over the extrapolations in the quark masses and those to the continuum limit.

## 4 The mass of the strange quark

Quark masses are important input parameters in many theoretical applications but, despite a great deal of activity, relatively little is known about their absolute values . Here I shall describe calculations of current quark masses defined through the PCAC relation. For charged kaons it can be written as

$$f_\mathrm{K}m_\mathrm{K}^2=(\overline{m}_\mathrm{u}+\overline{m}_\mathrm{s})0|\overline{u}\gamma _5s|\mathrm{K},$$
(10)
where the bars above the quark masses indicate that the “running” masses are considered, which depend on the energy scale. In order to determine the sum of up and strange quark masses using the experimental result for $`f_\mathrm{K}m_\mathrm{K}^2`$ one simply has to compute the matrix element $`0|\overline{u}\gamma _5s|\mathrm{K}`$ in a lattice simulation. However, in order to represent meaningful theoretical input, quark masses have to be renormalized, and hence their values depend on the adopted renormalization procedure. By convention quark masses are quoted in the $`\overline{\mathrm{MS}}`$ scheme of dimensional regularization at a reference scale, typically $`\mu =2\mathrm{GeV}`$. This implies that the renormalized matrix element in the $`\overline{\mathrm{MS}}`$ scheme has to be determined; it is related to its lattice counterpart by

$$0|\overline{u}\gamma _5s(\overline{\mu })|\mathrm{K}_{\overline{\mathrm{MS}}}=Z_\mathrm{P}^{\overline{\mathrm{MS}}}(g_0,a\overline{\mu })0|\overline{u}\gamma _5s|\mathrm{K}_{\mathrm{lat}}.$$
(11)
Here $`g_0`$ is the bare coupling, and $`\overline{\mu }`$ is the subtraction point in the $`\overline{\mathrm{MS}}`$ scheme. The renormalization factor $`Z_\mathrm{P}^{\overline{\mathrm{MS}}}`$ is known to one-loop order in lattice perturbation theory. However, it is well known that lattice perturbation theory converges slowly, and in order to remove all doubts about the reliability of the matching procedure it is evident that a non-perturbative determination of the renormalization factor is required.

However, a non-perturbative normalization condition relating the matrix elements in lattice regularization and the $`\overline{\mathrm{MS}}`$ scheme cannot be formulated, since the latter is only defined to any given order of perturbation theory. This technical difficulty can be overcome by introducing an intermediate scheme $`\mathrm{X}`$ and considering a two-step matching procedure, as shown schematically in Fig. 2. The first part is the matching of the bare current quark mass $`m_{\mathrm{lat}}(a)`$ to the running mass in the intermediate scheme $`\overline{m}_\mathrm{X}(\mu )`$. This amounts to computing $`Z_\mathrm{P}^\mathrm{X}(g_0,a\mu _0)`$ between the lattice scheme and scheme X at a fixed scale $`\mu _0`$. The second part is the determination of the scale dependence of the running mass $`\overline{m}_\mathrm{X}(\mu )`$ from $`\mu _0`$ up to very high energies, where the perturbative relation between $`\overline{m}_\mathrm{X}`$ and $`\overline{m}_{\overline{\mathrm{MS}}}`$ is expected to be reliable. Through this two-step process the use of lattice perturbation theory is completely avoided.

Two examples for such intermediate schemes have been proposed. The so-called “regularization-independent” (RI) scheme is described in . Another intermediate scheme is defined using the Schrödinger Functional (SF) of QCD . This scheme allows one to compute the scale dependence non-perturbatively over several orders of magnitude, using a recursive finite-size scaling technique. Once the scale dependence of the running mass in the SF scheme, $`\overline{m}_{\mathrm{SF}}`$, is known up to $`\mu 100\mathrm{GeV}`$ one can continue the scale evolution to infinite energy using the perturbative renormalization group functions and thereby extract the renormalization group invariant (RGI) quark mass $`M`$. At this point the matching to the $`\overline{\mathrm{MS}}`$ scheme is trivial, since $`M`$ is scheme-independent, i.e. it coincides in the SF and $`\overline{\mathrm{MS}}`$ schemes.

Figure 3 shows the non-perturbatively determined scale dependence of $`\overline{m}_{\mathrm{SF}}/M`$ computed in quenched QCD  and compares it to the perturbative scale evolution. The left-most data point in Fig. 3 corresponds to a scale $`\mu _0=275\mathrm{MeV}`$. It is here that the matching between lattice regularization and the $`\overline{\mathrm{MS}}`$ scheme via the SF is completed, through a non-perturbative calculation of $`Z_\mathrm{P}^{\mathrm{SF}}(g_0,a\mu _0)`$ at fixed $`\mu _0=275\mathrm{MeV}`$. All dependence on the intermediate SF scheme and the scale $`\mu _0`$ drops out in the total renormalization factor

$$\frac{\overline{m}_{\mathrm{SF}}(\mu _0)}{M}\times Z_\mathrm{P}^{\mathrm{SF}}(g_0,a\mu _0),$$
(12)
where $`\overline{m}_{\mathrm{SF}}(\mu _0)/M=0.864\pm 0.011`$, as can be read off the left-most point in Fig. 3 . The total renormalization factor is thus known with a precision of 1.5%.

### 4.1 The strange quark mass in quenched QCD

I shall now discuss the application of non-perturbative quark mass renormalization to compute the mass of the strange quark in the quenched approximation. Before discussing the details, it is useful to recall that mass ratios of light quarks are predicted quite accurately by Chiral Perturbation Theory (ChPT), for instance

$$\frac{M_\mathrm{s}}{\widehat{M}}=24.4\pm 1.4,\widehat{M}=\frac{1}{2}(M_\mathrm{u}+M_\mathrm{d}).$$
(13)
In order to determine $`M_\mathrm{s}`$, it is then sufficient to compute the sum $`(M_\mathrm{s}+\widehat{M})`$ on the lattice and combine it with eq. (13).

Here I shall concentrate on the calculation performed by the ALPHA/UKQCD Collaboration . By combining the PCAC relation, eq. (10), with the total renormalization factor, one obtains the sum of RGI quark masses $`(M_\mathrm{s}+\widehat{M})`$ in units of the kaon decay constant222Here, $`\mathrm{}`$ denotes the Dirac spinor for the isospin-symmetric light quark.

$$\frac{M_\mathrm{s}+\widehat{M}}{f_\mathrm{K}}=\left(\frac{M}{\overline{m}_{\mathrm{SF}}}\frac{1}{Z_\mathrm{P}}\right)\times \frac{m_\mathrm{K}^2}{0|\overline{\mathrm{}}\gamma _5s|\mathrm{K}_{\mathrm{lat}}}+O(a^2).$$
(14)
The dependence on the lattice spacing can be eliminated by performing a continuum extrapolation as shown in Fig. 4.

Using the experimental value $`f_\mathrm{K}=160\pm 2\mathrm{MeV}`$  one obtains the result for the sum of RGI quark masses in the continuum limit:

$$M_\mathrm{s}+\widehat{M}=140\pm 5\mathrm{MeV}.$$
(15)
This result can now be converted into $`\overline{m}_\mathrm{s}^{\overline{\mathrm{MS}}}(\overline{\mu })`$ at $`\overline{\mu }=2\mathrm{GeV}`$. First, eq. (15) is combined with the prediction from ChPT, eq. (13). Then, the relation between the running mass in the $`\overline{\mathrm{MS}}`$ scheme and the RGI quark mass is computed in 4-loop perturbation theory as

$$\overline{m}^{\overline{\mathrm{MS}}}(\overline{\mu })/M=0.7208\text{at}\overline{\mu }=2\mathrm{GeV}.$$
(16)
This procedure yields the final result in the quenched approximation

$$\overline{m}_\mathrm{s}^{\overline{\mathrm{MS}}}(2\mathrm{GeV})=97\pm 4\mathrm{MeV}.$$
(17)
The quoted uncertainty of $`\pm 4\mathrm{MeV}`$ contains all errors, except those due to quenching. The high precision of this result is a direct consequence of recent progress in lattice calculations, in particular the implementation of non-perturbative renormalization and control over lattice artefacts. As mentioned in section 2, the conversion into physical units is ambiguous in the quenched approximation. For $`\overline{m}_\mathrm{s}^{\overline{\mathrm{MS}}}(2\mathrm{GeV})`$ the resulting uncertainty was estimated in ref.  to amount to $`10\%`$.

Table 1 contains a compilation of recent results for $`m_\mathrm{s}`$ and $`\widehat{m}`$ obtained in the quenched approximation. Direct comparisons of these results should be made with care, since systematic errors have not been estimated in a uniform manner. Also, the conversion into physical units has been performed using different quantities. This manifests itself in the typical spread of the central values, which is of the order of 10%, consistent with the above estimate. It is still quite remarkable that lattice estimates for light quark masses have stabilized, which represents a big improvement with respect to the situation before 1998.

### 4.2 Sea quark effects in $`m_\mathrm{s}`$ and $`\widehat{m}`$

So far the most comprehensive study of quark masses using simulations with dynamical quarks has been presented recently by CP-PACS . Earlier results can be found in refs. . When discussing the results, it is important to keep in mind that non-perturbative renormalization has not yet been implemented for $`N_\mathrm{f}=2`$ flavours. Instead the renormalization of the lattice matrix element is usually performed using so-called “mean-field improved” lattice perturbation theory . In addition to computing the strange quark mass, CP-PACS have also estimated the average light quark mass through extrapolations to the chiral limit. Furthermore they have compared different definitions of the quark mass, based either on the axial vector Ward identity (AWI, i.e. the PCAC relation) or the vector Ward identity (VWI).

The continuum extrapolations for the strange quark mass are shown in Fig. 5. Results from the AWI and VWI definitions of the quark mass have been extrapolated enforcing a common continuum limit. By comparing the results for $`N_\mathrm{f}=2`$ to those obtained in the quenched approximation, one finds that the inclusion of dynamical quark effects decreases the estimate for the strange quark mass by roughly 20%. The results for $`m_\mathrm{s}`$ by CP-PACS can be summarized as follows:

$$m_\mathrm{s}^{\overline{\mathrm{MS}}}(2\mathrm{GeV})=\{\begin{array}{cc}\hfill 88_6^{+4}\mathrm{MeV},& N_\mathrm{f}=2\hfill \\ \hfill 110_4^{+3}\mathrm{MeV},& \text{quenched\hspace{0.17em}.}\hfill \end{array}$$
(18)
The light quark mass $`\widehat{m}`$ shows a similar decrease when sea quark effects are included. However, for both $`N_\mathrm{f}=2`$ and $`N_\mathrm{f}=0`$ (quenched) one finds that $`m_\mathrm{s}/\widehat{m}25`$ in the continuum limit, which is consistent with the prediction from ChPT, eq. (13).

Figure 5 also illustrates that lattice artefacts are quite substantial, so that an estimation of quark masses at non-zero lattice spacing would be misleading. This underlines once more the importance of the continuum extrapolation.

## 5 Concluding remarks

Lattice simulations of quenched QCD have reached a level of precision of a few per cent, thanks to increases in computer power and conceptual advances. The picture that has emerged is that the quenched approximation works surprisingly well, as signified by the success in predicting the spectrum of low-lying hadrons. Non-perturbative renormalization and good control over lattice artefacts are the crucial ingredients that lead to precise estimates of the mass of the strange quark. Also, in order to control the region of very light quarks, which cannot be simulated on present machines, it is useful to combine lattice simulations with Chiral Perturbation Theory. In fact, a combination of these techniques may be instrumental in resolving the long-standing question of whether the up-quark is massless . A pilot study along these lines has already appeared in the literature .

A major share of computational resources is devoted to studying the effects of dynamical quarks. First results indicate that sea quarks “are doing the right thing”, since the discrepancies in the light hadron spectrum between simulation and experiment are decreased. However, more effort is required to understand sea quark effects on a quantitative level. This not only requires larger computers but also the implementation of additional technology such as non-perturbative renormalization for the dynamical case. Finally, there are algorithmic challenges: the efficient simulation of sea quarks for which $`m_{\mathrm{PS}}/m_\mathrm{V}0.5`$, and also of odd $`N_\mathrm{f}`$ is crucial for the development of lattice QCD into a fully fledged phenomenological tool.

I am grateful to the organizers, in particular Dr Grazyna Odyniec, for the invitation and the offered support. I wish to thank all organizers for creating such a stimulating and pleasant environment at the conference.

## References