text:
# Bogoliubov Approach to Superfluidity of Atoms in an Optical Lattice

## 1 Introduction

Spectacular progress has been made in experimental studies of atoms loaded into an optical lattice in the region of the Mott superfluid insulator quantum phase transition . In this article, we shall discuss the superfluid density and the quasi-momentum distribution, which is directly related to the matter-wave interference patterns that can be observed in such experiments. To do this we use the Bogoliubov method , as developed for use in optical lattices . In a previous paper  we used this method to produce results for squeezing that are consistent with those of other approaches previously reported in the literature . In this paper we shall show how it can be used to predict the decrease in the superfluid fraction and the corresponding variations in the matter wave interference fringes that should be directly observable in future experiments. This extends our previous studies based on exact calculation for small one-dimensional systems  into the experimentally relevant regime of lattice sizes and particle numbers.

We first introduce the Bose-Hubbard Hamiltonian for atoms in an optical lattice . We then describe briefly our method for determining the superfluid fraction based on the rigidity of the system under a twist of the condensate phase . Using a perturbative formulation analogous to the Drude weight  the Bogoliubov approximation gives us a particularly direct way of finding this quantity. It also gives a simple picture of how superfluidity is suppressed by quantum depletion of the condensate. We shall compare the results for various quantities calculated using the Bogoliubov approximation with exact numerical calculations for the case of modest numbers of atoms and lattice sites .

## 2 The Bose-Hubbard Model and Superfluidity

The Bose-Hubbard Hamiltonian for atoms in an one-dimensional optical lattice with $`I`$ sites has the form :

$$\widehat{H}=\underset{i=1}{\overset{I}{}}\widehat{n}_iϵ_iJ\underset{i=1}{\overset{I}{}}(\widehat{a}_{i+1}^{}\widehat{a}_i+\widehat{a}_i^{}\widehat{a}_{i+1})+\frac{V}{2}\underset{i=1}{\overset{I}{}}\widehat{n}_i(\widehat{n}_i1).$$
(1)
Here $`J`$ represents the coupling between adjacent lattice sites due to tunneling and $`V`$ is the strength of repulsion between atoms on the same site. The non-interacting energy of the atoms on each site, $`ϵ_i`$, will have some variation that is typically smooth on the scale of the condensate. We shall consider below both, the case where this is a constant, as well as the extension to the case where it varies in a trapped condensate. This Bose-Hubbard Hamiltonian should be an appropriate model when the loading process produces atoms in the lowest vibrational state of each well, with a chemical potential smaller than the distance to the first vibrationally excited state. This is known to be possible from the results of recent experiments .

The concept of superfluidity is closely related to the existence of a condensate in the interacting many–body system. Formally, the one–body density matrix $`\rho ^{(1)}(\stackrel{}{x},\stackrel{}{x}^{})`$ has to have exactly one macroscopic eigenvalue which defines the number of particles in the condensate; the corresponding eigenvector describes the condensate wave function $`\varphi _0\left(\stackrel{}{x}\right)=e^{i\mathrm{\Theta }(\stackrel{}{x})}\left|\varphi _0\left(\stackrel{}{x}\right)\right|`$. A spatially varying condensate phase, $`\mathrm{\Theta }\left(\stackrel{}{x}\right)`$, is associated with a velocity field for the condensate by

$$\stackrel{}{v}_0\left(\stackrel{}{x}\right)=\frac{\mathrm{}}{m}\stackrel{}{}\mathrm{\Theta }\left(\stackrel{}{x}\right).$$
(2)
This irrotational velocity field is identified with the velocity of the *superfluid flow*, $`\stackrel{}{v}_s\left(\stackrel{}{x}\right)\stackrel{}{v}_0\left(\stackrel{}{x}\right)`$ (,) and enables us to derive an expression for the superfluid fraction, $`f_s`$. Consider a system with a finite linear dimension, $`L`$, in the $`\stackrel{}{e}_1`$–direction and a ground–state energy, $`E_0`$, calculated with periodic boundary conditions. Now we impose a linear phase variation, $`\mathrm{\Theta }\left(\stackrel{}{x}\right)=\theta x_1/L`$ with a total twist angle $`\theta `$ over the length of the system in the $`\stackrel{}{e}_1`$–direction. The resulting ground–state energy, $`E_\theta `$ will depend on the phase twist. For very small twist angles, $`\theta \pi `$, the energy difference, $`E_\theta E_0`$, can be attributed to the kinetic energy, $`T_s`$, of the superflow generated by the phase gradient. Thus,

$$E_\theta E_0=T_s=\frac{1}{2}mNf_s\stackrel{}{v}_s^2,$$
(3)
where $`m`$ is the mass of a single particle and $`N`$ is the total number of particles so that $`mNf_s`$ is the total mass of the superfluid component. Replacing the superfluid velocity, $`\stackrel{}{v}_s`$ with the phase gradient according to Eq. (2) leads to a fundamental relation for the superfluid fraction

$$f_s=\frac{2m}{\mathrm{}^2}\frac{L^2}{N}\frac{E_\theta E_0}{\theta ^2}=\frac{1}{N}\frac{E_\theta E_0}{J\left(\mathrm{\Delta }\theta \right)^2},$$
(4)
where the second equality applies to a lattice system on which a linear phase variation has been imposed. Here the distance between sites is $`a`$, the phase variation over this distance is $`\mathrm{\Delta }\theta `$, and the number of sites is $`I`$. In this case, $`J\mathrm{}^2/(2ma^2)`$.

Technically the phase variation can be imposed through so-called twisted boundary conditions . In the context of the discrete Bose-Hubbard model it is, however, more convenient to map the phase variation by means of a unitary transformation onto the Hamiltonian. The resulting “twisted” Hamiltonian

$$\widehat{H}_\theta =\underset{i=1}{\overset{I}{}}\widehat{n}_iϵ_iJ\underset{i=1}{\overset{I}{}}(e^{\mathrm{i}\mathrm{\Delta }\theta }\widehat{a}_{i+1}^{}\widehat{a}_i+e^{\mathrm{i}\mathrm{\Delta }\theta }\widehat{a}_i^{}\widehat{a}_{i+1})+\frac{V}{2}\underset{i=1}{\overset{I}{}}\widehat{n}_i(\widehat{n}_i1)$$
(5)
exhibits additional phase factors $`e^{\pm \mathrm{i}\mathrm{\Delta }\theta }`$ — the so-called Peierls phase factors — in the hopping term . These phase factors show that the twist is equivalent to the imposition of an acceleration on the lattice for a finite time. It is interesting to note that the present experiments enable us to make a specific connection between the formal and operational aspects of the system.

We calculate the change in energy $`E_\theta E_0`$ under the assumption that the phase change $`\mathrm{\Delta }\theta `$ is small so that we can write:

$$e^{\mathrm{i}\mathrm{\Delta }\theta }1\mathrm{i}\mathrm{\Delta }\theta \frac{1}{2}(\mathrm{\Delta }\theta )^2.$$
(6)
Using this expansion the twisted Hamiltonian (5) takes the following form:

$$\widehat{H}_\theta \widehat{H}_0+\mathrm{\Delta }\theta \widehat{J}\frac{1}{2}(\mathrm{\Delta }\theta )^2\widehat{T}=\widehat{H}_0+\widehat{H}_{\text{pert}},$$
(7)
where we retain terms up to second order in $`\mathrm{\Delta }\theta `$. The current operator $`\widehat{J}`$ (N.B. that the physical current is given by this expression multiplied by $`\frac{1}{\mathrm{}}`$) and the hopping operator $`\widehat{T}`$ are given by:

$`\widehat{J}`$ $`=`$ $`\mathrm{i}J{\displaystyle \underset{i=1}{\overset{I}{}}}(\widehat{a}_{i+1}^{}\widehat{a}_i\widehat{a}_i^{}\widehat{a}_{i+1})`$ (8)
$`\widehat{T}`$ $`=`$ $`J{\displaystyle \underset{i=1}{\overset{I}{}}}(\widehat{a}_{i+1}^{}\widehat{a}_i+\widehat{a}_i^{}\widehat{a}_{i+1}).`$ (9)
The change in the energy $`E_\theta E_0`$ due to the imposed phase twist can now be evaluated in second order perturbation theory

$$E_\theta E_0=\mathrm{\Delta }E^{(1)}+\mathrm{\Delta }E^{(2)}.$$
(10)
The first order contribution to the energy change is proportional to the expectation value of the hopping operator

$$\mathrm{\Delta }E^{(1)}=\mathrm{\Psi }_0|\widehat{H}_{\text{pert}}|\mathrm{\Psi }_0=\frac{1}{2}(\mathrm{\Delta }\theta )^2\mathrm{\Psi }_0|\widehat{T}|\mathrm{\Psi }_0.$$
(11)
Here $`|\mathrm{\Psi }_0`$ is the ground state of the original Bose-Hubbard Hamiltonian (1). The second order term is related to the matrix elements of the current operator involving the excited states $`|\mathrm{\Psi }_\nu `$ ($`\nu =1,2,\mathrm{}`$) of the original Hamiltonian

$$\mathrm{\Delta }E^{(2)}=\underset{\nu 0}{}\frac{|\mathrm{\Psi }_\nu |\widehat{H}_{\text{pert}}|\mathrm{\Psi }_0|^2}{E_\nu E_0}=(\mathrm{\Delta }\theta )^2\underset{\nu 0}{}\frac{|\mathrm{\Psi }_\nu |\widehat{J}|\mathrm{\Psi }_0|^2}{E_\nu E_0}.$$
(12)
Thus we obtain for the energy change up to second order in $`\mathrm{\Delta }\theta `$

$`E_\theta E_0`$ $`=`$ $`(\mathrm{\Delta }\theta )^2\left({\displaystyle \frac{1}{2}}\mathrm{\Psi }_0|\widehat{T}|\mathrm{\Psi }_0{\displaystyle \underset{\nu 0}{}}{\displaystyle \frac{|\mathrm{\Psi }_\nu |\widehat{J}|\mathrm{\Psi }_0|^2}{E_\nu E_0}}\right)=I(\mathrm{\Delta }\theta )^2D,`$
$`D`$ $``$ $`{\displaystyle \frac{1}{I}}\left({\displaystyle \frac{1}{2}}\mathrm{\Psi }_0|\widehat{T}|\mathrm{\Psi }_0{\displaystyle \underset{\nu 0}{}}{\displaystyle \frac{|\mathrm{\Psi }_\nu |\widehat{J}|\mathrm{\Psi }_0|^2}{E_\nu E_0}}\right).`$ (13)
The quantity $`D`$, defined above, is formally equivalent to the Drude weight used to specify the DC conductivity of charged fermionic systems . The superfluid fraction is then given by the contribution of both the first and second order term:

$`f_\mathrm{s}`$ $`=`$ $`f_\mathrm{s}^{(1)}f_\mathrm{s}^{(2)};`$
$`f_\mathrm{s}^{(1)}`$ $``$ $`{\displaystyle \frac{1}{2NJ}}\left(\mathrm{\Psi }_0|\widehat{T}|\mathrm{\Psi }_0\right),`$ (14)
$`f_\mathrm{s}^{(2)}`$ $``$ $`{\displaystyle \frac{1}{NJ}}\left({\displaystyle \underset{\nu 0}{}}{\displaystyle \frac{|\mathrm{\Psi }_\nu |\widehat{J}|\mathrm{\Psi }_0|^2}{E_\nu E_0}}\right).`$
Here $`N`$ is the number of atoms in the lattice. In general both, the first and the second order term contribute. For a translationally invariant lattice the second term vanishes (as is going to be shown latter) in the Bogoliubov limit that we shall use in this study. However, in exact calculations and in the Bogoliubov approximation for an inhomogeneous lattice the second order term plays a role.

We can further understand this approach to the superfluid density by calculating the flow that is produced by the application of the phase twist. To do this we work out the expectation value of the current operator expressed in terms of the twisted variables:

$$\widehat{J}_\theta =\mathrm{i}J\underset{i=1}{\overset{I}{}}(e^{\mathrm{i}\mathrm{\Delta }\theta }\widehat{a}_{i+1}^{}\widehat{a}_ie^{\mathrm{i}\mathrm{\Delta }\theta }\widehat{a}_i^{}\widehat{a}_{i+1}).$$
(15)
We expand this to find the lowest order contributions, i.e.:

$$\widehat{J}_\theta \widehat{J}+J\mathrm{\Delta }\theta \underset{i=1}{\overset{I}{}}(\widehat{a}_{i+1}^{}\widehat{a}_i+\widehat{a}_i^{}\widehat{a}_{i+1})=\widehat{J}\widehat{T}\mathrm{\Delta }\theta .$$
(16)
We use first order perturbation theory on the wave function to obtain the following expression:

$`\mathrm{\Psi }(\mathrm{\Delta }\theta )|\widehat{J}_\theta |\mathrm{\Psi }(\mathrm{\Delta }\theta )`$ $`=`$ $`2\mathrm{\Delta }\theta \left({\displaystyle \frac{1}{2}}\mathrm{\Psi }_0|\widehat{T}|\mathrm{\Psi }_0{\displaystyle \underset{\nu 0}{}}{\displaystyle \frac{|\mathrm{\Psi }_\nu |\widehat{J}|\mathrm{\Psi }_0|^2}{E_\nu E_0}}\right)`$ (17)
$`=`$ $`2NJf_\mathrm{s}\mathrm{\Delta }\theta .`$ (18)
If we note that the kinetic energy for a small quasi-momentum $`q`$ on a lattice is given by $`Jq^2a^2`$, we can define the effective mass as $`m^{}=\frac{\mathrm{}^2}{2Ja^2}`$. Here the quasi-momenta are given by $`q=\frac{2\pi }{Ia}j`$ with $`j=1,\mathrm{},(I1)`$ and lattice spacing $`a`$. Thus, the physical current, Eq. (18) multiplied by $`\frac{1}{\mathrm{}}`$, can be expressed as:

$$\mathrm{\Psi }(\mathrm{\Delta }\theta )|\widehat{J}_\theta |\mathrm{\Psi }(\mathrm{\Delta }\theta )=Nf_\mathrm{s}\mathrm{\Delta }\theta \frac{\mathrm{}}{m^{}a^2}.$$
(19)
This is the total flux and we need to divide $`I`$ to get the flux density, i.e.

$`{\displaystyle \frac{1}{I}}\mathrm{\Psi }(\mathrm{\Delta }\theta )\left|\widehat{J}_\theta \right|\mathrm{\Psi }(\mathrm{\Delta }\theta )`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{}\mathrm{\Delta }\theta }{m^{}a}}\right)\left({\displaystyle \frac{Nf_\mathrm{s}}{aI}}\right)`$ (20)
$`=`$ $`v_\text{s}n_\text{s}.`$
So we see that the Drude formulation of the superfluid fraction (2) gives an intuitively satisfying expression for the amount of flowing superfluid.

## 3 The Bogoliubov Approximation to the Bose-Hubbard Hamiltonian

We use the Bogoliubov approximation for the Bose-Hubbard model in the limit that quantum fluctuations, or equivalently depletion of the condensate, is not too great. In the limit that the quantum depletion can be completely ignored, we can replace the creation and annihilation operators, $`\widehat{a}_i^{}`$ and $`\widehat{a}_i`$, on each site with a $`c`$-number, $`z_i`$. This leads to a set of coupled non-linear Schrödinger, i.e. Gross-Pitaevskii (GP), equations for these amplitudes :

$$\mathrm{i}\mathrm{}_tz_i=J(z_{i+1}+z_{i1})+Vz_iz_i^{}z_i.$$
(21)
This equation can be used to study the properties of the condensate loaded into the lattice when the tunneling kinetic energy is large enough compared to the interaction energy though small enough for the one-band Bose-Hubbard model to be valid. We then include the quantum fluctuations in our description of the system using the Bogoliubov approximation, where we suppose that we can write the full annihilation operator in terms of the $`c`$-number part and a fluctuation operator thus:

$$\widehat{a}_i=(z_i+\widehat{\delta }_i)e^{\mathrm{i}\frac{\mu t}{\mathrm{}}}.$$
(22)
This form will be useful when we are looking at the properties of a time-independent or adiabatic ground state. In using this method we are assuming that the fluctuation part is small. The Bogoliubov method gives us expressions for the averages of the squares of the fluctuation operator and allows us to determine whether this assumption is valid. We shall examine its validity by comparing the results for various physical quantities with exact numerical calculations based on the Bose-Hubbard Hamiltonian.

### 3.1 Bogoliubov theory for the translationally invariant lattice

The ground state solution for the translationally invariant lattice gives the eigenvalue:

$$\mu =n_0V2J,$$
(23)
where

$$n_0=N/I$$
(24)
is the mean number of atoms on each site of the lattice. We take $`N`$ to be the total number of atoms and $`I`$ to be the number of sites in the one-dimensional lattice.

The Bogoliubov equations for the lattice have the following form:

$$\mathrm{i}\mathrm{}_t\widehat{\delta }_i=(2n_0V\mu )\widehat{\delta }_iJ(\widehat{\delta }_{i+1}+\widehat{\delta }_{i1})+n_0V\widehat{\delta }_i^{}.$$
(25)
This is solved by constructing quasi-particles for the lattice which diagonalize the Hamiltonian , i.e

$`\widehat{\delta }_i`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{I}}}{\displaystyle \underset{q}{}}[u^q\widehat{\alpha }_qe^{\mathrm{i}(qia\omega _qt)}v^q\widehat{\alpha }_q^{}e^{\mathrm{i}(qia\omega _qt)}]`$ (26)
$`\widehat{\delta }_i^{}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{I}}}{\displaystyle \underset{q}{}}[u^q\widehat{\alpha }_q^{}e^{\mathrm{i}(qia\omega _qt)}v^q\widehat{\alpha }_qe^{\mathrm{i}(qia\omega _qt)}],`$ (27)
where $`a`$ is the lattice spacing. The quasi-particle operators obey the usual Bose commutation relations:

$$[\widehat{\alpha }_q,\widehat{\alpha }_q^{}^{}]=\delta _{qq^{}}$$
(28)
and have the following expectation values at some temperature $`T`$:

$$\widehat{\alpha }_q^{}\widehat{\alpha }_q^{}=\delta _{qq^{}}[\mathrm{exp}(\mathrm{}\omega _q/k_bT)1]^1.$$
(29)
We then find the following equations for the excitation amplitudes and frequencies,

$`\mathrm{}\omega _qu^q`$ $`=`$ $`\left[n_0V+4J\mathrm{sin}^2\left({\displaystyle \frac{qa}{2}}\right)\right]u^qn_0Vv^q,`$ (30)
$`\mathrm{}\omega _qv^q`$ $`=`$ $`\left[n_0V+4J\mathrm{sin}^2\left({\displaystyle \frac{qa}{2}}\right)\right]v^qn_0Vu^q.`$ (31)
Thus, the expressions for the $`u^q`$ and $`v^q`$ yield:

$`|u^q|^2`$ $`=`$ $`{\displaystyle \frac{K(q)+n_0V+\mathrm{}\omega _q}{2\mathrm{}\omega _q}}`$ (32)
$`|v^q|^2`$ $`=`$ $`{\displaystyle \frac{K(q)+n_0V\mathrm{}\omega _q}{2\mathrm{}\omega _q}},`$ (33)
where the phonon excitation frequencies are given by:

$`\mathrm{}\omega _q`$ $`=`$ $`\sqrt{K(q)[2n_0V+K(q)]}`$ (34)
$`K(q)`$ $`=`$ $`4J\mathrm{sin}^2\left({\displaystyle \frac{qa}{2}}\right).`$ (35)
### 3.2 Expressions for the number superfluid fraction in the translationally invariant lattice

Having obtained the expressions for the excitations we can now determine the superfluid fraction. The quantity we need to calculate is just the first order term of the Drude expression (2), because the second order term vanishes in the Bogoliubov limit due to the translational invariance of the lattice \[see Eq. (48)\], i.e.

$$f_\mathrm{s}=\frac{1}{2NJ}\mathrm{\Psi }_0|\widehat{T}|\mathrm{\Psi }_0=\frac{1}{2N}\underset{i=1}{\overset{I}{}}\mathrm{\Psi }_0|\widehat{a}_{i+1}^{}\widehat{a}_i+\widehat{a}_i^{}\widehat{a}_{i+1}|\mathrm{\Psi }_0.$$
(36)
In the Bogoliubov approximation this has the form:

$`f_\mathrm{s}`$ $`=`$ $`{\displaystyle \frac{1}{2N}}{\displaystyle \underset{i=1}{\overset{I}{}}}\mathrm{\Psi }_0|(\widehat{\delta }_{i+1}^{}+z_{i+1})(\widehat{\delta }_i+z_i)+(\widehat{\delta }_i^{}+z_i)(\widehat{\delta }_{i+1}+z_{i+1})|\mathrm{\Psi }_0`$ (37)
$`=`$ $`{\displaystyle \frac{1}{2N}}{\displaystyle \underset{i=1}{\overset{I}{}}}\mathrm{\Psi }_0|2z_i^2+\widehat{\delta }_{i+1}^{}\widehat{\delta }_i+\widehat{\delta }_i^{}\widehat{\delta }_{i+1}|\mathrm{\Psi }_0.`$
We can now express the fluctuation operators, Eqs. (26) and (27), in terms of the quasi-particle operators that diagonalize the quadratic Hamiltonian. This leads to the following expression for the superfluid fraction at finite temperature:

$`f_\mathrm{s}`$ $`=`$ $`{\displaystyle \frac{1}{2N}}[{\displaystyle \underset{i=1}{\overset{I}{}}}2z_i^2+{\displaystyle \frac{1}{I}}{\displaystyle \underset{q}{}}[u^q\widehat{\alpha }_qe^{\mathrm{i}q(i+1)a}v^q\widehat{\alpha }_q^{}e^{\mathrm{i}q(i+1)a}]`$ (38)
$`\times {\displaystyle \underset{q^{}}{}}[u^q^{}\widehat{\alpha }_q^{}^{}e^{\mathrm{i}q^{}ia}v^q^{}\widehat{\alpha }_q^{}e^{\mathrm{i}q^{}ia}]`$
$`+{\displaystyle \frac{1}{I}}{\displaystyle \underset{i=1}{\overset{I}{}}}{\displaystyle \underset{q}{}}[u^q\widehat{\alpha }_q^{}e^{\mathrm{i}qia}v^q\widehat{\alpha }_qe^{\mathrm{i}qia}]`$
$`\times {\displaystyle \underset{q^{}}{}}[u^q^{}\widehat{\alpha }_q^{}e^{\mathrm{i}q^{}(i+1)a}v^q^{}\widehat{\alpha }_q^{}^{}e^{\mathrm{i}q^{}(i+1)a}]],`$
and we find in the zero temperature limit of a translationally invariant lattice:

$$f_\mathrm{s}=\frac{I}{N}\left[z^2+\frac{1}{I}\underset{q}{}|v^q|^2\mathrm{cos}(qa)\right].$$
(39)
Here the summation runs over all quasi-momenta $`q=\frac{2\pi }{Ia}j`$ with $`j=1,\mathrm{},(I1)`$ and we have called $`z`$ the value of all $`z_i`$ in a translationally invariant system. This shows that in the limit of zero lattice spacing (while keeping $`q`$ finite) the superfluid fraction is unity as we have the normalization condition:

$$Iz^2+\underset{q}{}|v^q|^2=N.$$
(40)
These expressions give a direct insight into the change of the superfluid fraction as atoms are pushed out of the condensate due to interactions. In Eq. (39) the sum involving the Bogoliubov amplitudes $`v^q`$ characterizes the difference between the condensate fraction, which is given by the first term, and the superfluid fraction. For weak interactions and a small depletion, which fills only the lower quarter of the band where the $`\mathrm{cos}(qa)`$ term has a positive sign, the superfluid fraction is larger than the condensate fraction. Thus the depletion of the condensate has initially little effect on superfluidity. When the depleted population spreads into the central part of the band, where the $`\mathrm{cos}(qa)`$ term has a negative sign, the superfluid fraction is reduced and might even become smaller than the condensate fraction. Finally, the population in the upper quarter of the band again produces a positive contribution to the superflow. In a sense the interactions are playing a role akin to Fermi exclusion ”pressure” in the case of electron flow in a band. This, however can lead to perfect filling and cancellation of the flow. In the case of our Bogoliubov description we can only see reduction of the flow, not a perfect switching off of the superfluid. This happens in the Mott insulator state, which cannot be described by the Bogoliubov approximation.

In the next section we outline the version of the Bogoliubov theory that should be best suited to treating these problems, i.e self-consistent Bogoliubov theory.

## 4 Self-consistent HFB-Popov theory

In this section we explore the limits of validity of the simplest zero temperature self-consistent Bogoliubov theory, a simplified version of the Hartree-Fock-Bogoliubov approximation originally introduced by Popov . The HFB-Popov theory is an extension of the standard Bogoliubov approximation that takes into account the depletion of the condensate but neglects the anomalous average. As discussed in the previous section, taking into account the depletion of the condensate is important as the transition is approached because the depleted population causes the reduction of the superfluidity. Although the HFB-Popov aproach has the limitation that it doesn’t take into account the full effect of the medium because it neglects the anomalous average , it can be considered a better theory for the elementary excitations than the full HFB due to the fact that the theory is gapless and doesn’t violate Goldstone’s theorem.

A derivation of the Bogoliubov equations for the quasiparticle amplitudes in a lattice can be found for example in Ref. . These equations only take into account terms up to second order in the fluctuations. Including third and fourth order terms by treating them in a self-consistent mean field approximation ( ,  ) and neglecting anomalous average terms yields the following HFB-Popov equations:

$`\mathrm{}\omega _qu_i^q+c^qz_i`$ $`=`$ $`(2V(|z_i|^2+\stackrel{~}{n}_i)\mu +ϵ_i)u_i^qJ(u_{i+1}^q+u_{i1}^q)Vz_i^2v_i^q,`$ (41)
$`\mathrm{}\omega _qv_i^qc^qz_i^{}`$ $`=`$ $`(2V(|z_i|^2+\stackrel{~}{n}_i)\mu +ϵ_i)v_i^qJ(v_{i+1}^q+v_{i1}^q)Vz_i^2u_i^q,`$ (42)
$`\mu z_i`$ $`=`$ $`J(z_{i+1}+z_{i1})+\left(V(|z_i|^2+2\stackrel{~}{n}_i)+ϵ_i\right)z_i,`$ (43)
$`\stackrel{~}{n}_i`$ $`=`$ $`{\displaystyle \underset{q}{}}\left|v_i^q\right|^2,`$ (44)
$`N`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{I}{}}}(|z_i|^2+\stackrel{~}{n}_i),`$ (45)
$`c^q`$ $`=`$ $`V{\displaystyle \underset{i}{}}|z_i|^2(z_i^{}u_i^qz_iv_i^q).`$ (46)
Where $`ϵ_i`$ is the energy offset at site $`i`$ due to an external potential, $`\{u_i^q,v_i^q\}`$ and $`\omega _q`$ are respectively the quasiparticle amplitudes and energies, thus

$$\widehat{\delta }_i=\underset{q}{}u_i^q\widehat{\alpha }_qe^{\mathrm{i}\omega _qt}v_i^q\widehat{\alpha }_q^{}e^{\mathrm{i}\omega _qt},$$
(47)
$`\stackrel{~}{n}_i`$ is the depletion at site $`i`$, and $`N`$ is the total number of particles. The parameters $`c^q`$ ensure the $`\{u_i^q,v_i^q\}`$ solutions to the above equations with $`\omega _q0`$ to be orthogonal to the condensate (Ref. ).

By calculating the quasiparticle amplitudes and the condensate density it is possible to get information about most of the physical properties of the system. For example the superfluid fraction and the on site number fluctuations are given by:

$`f_\mathrm{s}`$ $`=`$ $`f_s^{(1)}f_s^{(2)},`$ (48)
$`f_s^{(1)}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{I}{}}}f_{si}^{(1)}={\displaystyle \frac{1}{2N}}{\displaystyle \underset{i=1}{\overset{I}{}}}\left[(z_{i+1}z_i^{}+z_{i+1}^{}z_i)+{\displaystyle \underset{q}{}}(v_i^qv_{i+1}^q+v_i^qv_{i+1}^q)\right],`$
$`f_s^{(2)}`$ $`=`$ $`{\displaystyle \frac{J}{N}}\left({\displaystyle \underset{q,q^{}}{}}{\displaystyle \frac{\left|_i\left(u_{i+1}^qv_i^q^{}u_i^qv_{i+1}^q^{}\right)\right|^2}{\mathrm{}\omega _q+\mathrm{}\omega _q^{}}}+\delta _{qq^{}}{\displaystyle \frac{\left|_i(u_{i+1}^qv_i^qu_i^qv_{i+1}^q)\right|^2}{2\mathrm{}\omega _q}}\right),`$
$$\mathrm{\Delta }n_i^2=|z_i|^2\underset{q}{}|u_i^qv_i^q|^2.$$
(49)
From the complete expression of the superfluid fraction, it can be seen explicitly how due to the translational invariance, the second order term vanishes in the homogeneous system.

### 4.1 Translationally Invariant lattice

For the translational invariant lattice we use the quasiparticle transformation given by Eqs. (26) and (27). Under this transformation the self consistent equations can be written, generalizing the previous version, as:

$`\mu `$ $`=`$ $`\left(|z|^2+{\displaystyle \frac{2}{I}}{\displaystyle \underset{q}{}}|v^q|^2\right)V2J,`$ (50)
$`|u^q|^2`$ $`=`$ $`{\displaystyle \frac{K(q)+|z|^2V+\mathrm{}\omega _q}{2\mathrm{}\omega _q}},`$ (51)
$`|v^q|^2`$ $`=`$ $`{\displaystyle \frac{K(q)+|z|^2V\mathrm{}\omega _q}{2\mathrm{}\omega _q}},`$ (52)
$`N`$ $`=`$ $`I|z|^2+{\displaystyle \underset{q}{}}|v^q|^2.`$ (53)
Here the phonon excitation spectrum is given by:

$$\mathrm{}\omega _q=\sqrt{K(q)[2|z|^2V+K(q)]},$$
(54)
and $`K(q)`$ is given by Eq. (35). Again we omit the subscript in the amplitudes $`z_i`$ because they have the same value at all lattice sites. Notice also that due to translational invariance, the $`c^q`$ coefficients vanish.

In the homogeneous system, the form of the HFB-Popov equations for the quasiparticle amplitudes and energies is very close to the standard Bogoliubov form. We do, however, have to replace $`n_0=N/I`$ by the condensate amplitude $`|z|^2`$ which must take into account the depletion of the condensate. We solve for the condensate amplitude as a function of the external parameters $`J,V,N`$ and $`I`$ by inserting Eq. (52) in Eq. (53). Once $`|z|^2`$ is determined, we use it to calculate the other expressions.

In Fig. 1 we compare the number fluctuations on a lattice site, the condensate fraction and the total and second order superfluid fraction determined from the exact solution of the Bose-Hubbard Hamiltonian to the self consistent HFB-Popov predictions as a function of the ratio $`V_{eff}=V/J`$. The systems used for the comparisons have three wells, $`I=3`$, and commensurate filling factors $`n_0=5,10,20`$ and 50. Similar results for the incommensurate case with $`N=16,31,61,151`$ are shown in Fig. 2. We were restricted to consider only three wells due to computational limitations. The size of the matrix needed in the exact solution for $`N`$ atoms and $`I`$ wells scales as $`\frac{(N+I1)!}{N!(I1)!}`$. However, if the HFB-Popov approach works well for these small systems we expect it to provide a good description of the larger systems prepared in the lab.

Because the second order term of the superfluid fraction (second term of Eq.2) vanishes in the HFB-Popov approach (see Eq. 48), we only expect the self consistent HFB-Popov theory to give a good description of the superfluid fraction in the region where the second order term is extremely small, provided it predicts accurately the first order term. This is exactly what is observed in the plots. When the second order term starts to grow, typically above $`0.5V_{eff}^{_{\text{crit}}}`$, the HFB-Popov theory starts to fail. An estimate of $`V_{eff}^{\text{crit}}`$ is shown by a vertical line in some of the figures. This was obtained by using the second order perturbative approach presented in Ref. . With increasing filling factor the critical value is shifted towards larger values of the interaction strength, and the region in which the HFB-Popov theory is accurate gets larger. It is interesting to note that the number fluctuations predicted by the theory are accurate in a greater range than the other physical quantities shown. Its predictions of squeezing agree very well with the exact solutions right up to the point where the number fluctuations become less than unity.

For the cases with non-commensurate fillings depicted in Fig. 2, the agreement is significantly better for all quantities. This is not surprising because when the filling is not commensurate there is always a superfluid present and the Mott transition doesn’t occur. As can be seen in the plots for these cases the second order term is always very small.

### 4.2 Inhomogeneous lattice

In this section we consider the experimentally relevant case when there is an external magnetic confinement in addition to the lattice potential. In this situation, the self consistent HFB-Popov equations take the form:

$`\mathrm{}\omega _qu_i^q+c^qz_i`$ $`=`$ $`\left(2V\left(\left|z_i\right|^2+\stackrel{~}{n}_i\right)\mu +\mathrm{\Omega }i^2\right)u_i^qJ\left(u_{i+1}^q+u_{i1}^q\right)Vz_i^2v_i^q,`$ (55)
$`\mathrm{}\omega _qv_i^qc^qz_i^{}`$ $`=`$ $`\left(2V\left(\left|z_i\right|^2+\stackrel{~}{n}_i\right)\mu +\mathrm{\Omega }i^2\right)v_i^qJ\left(v_{i+1}^q+v_{i1}^q\right)Vz_i^2u_i^q,`$ (56)
$`\mu z_i`$ $`=`$ $`J\left(z_{i+1}+z_{i1}\right)+\left(V(|z_i|^2+2\stackrel{~}{n}_i)+\mathrm{\Omega }i^2\right)z_i,`$ (57)
$`\stackrel{~}{n}_i`$ $`=`$ $`{\displaystyle \underset{q}{}}|v_i^q|^2,`$ (58)
$`N`$ $`=`$ $`{\displaystyle \underset{i}{}}(|z_i|^2+\stackrel{~}{n}_i).`$ (59)
$`c^q`$ $`=`$ $`V{\displaystyle \underset{i}{}}|z_i|^2(z_i^{}u_i^qz_iv_i^q)`$ (60)
where $`\mathrm{\Omega }=\frac{1}{2}m\omega ^2a^2`$, with $`m`$ the mass of the atoms, $`\omega `$ the trap frequency, and $`a`$ the lattice spacing. The site indices $`i`$ are chosen such that $`i=0`$ corresponds to the center of the trap. Again the $`c^q`$are introduced to ensure the orthogonality of the excitations to the condensate . We have solved the HFB-Popov equations for this system by an iterative procedure, similar to the one followed in Ref. . Each cycle of the iteration consists of two steps. In the first step we solve Eq. (57) subject to the constraint Eq. (59) by using the $`\stackrel{~}{n}_i`$ obtained in the previous cycle. This generates new values for the $`z_i`$. In the second step we solve for $`\{u_i^q,v_i^q\}`$ in Eqs. (55) using the $`\stackrel{~}{n}_i`$ from the previous cycle and the newly generated $`z_i`$. The $`\{u_i^q,v_i^q\}`$ are used then to update $`\stackrel{~}{n}_i`$. Because the HFB-Popov is gapless, it is possible to keep the orthogonality of the excitations to the condensate by solving Eqs. (55) with the $`c^q`$ set to zero but removing in each cycle the projection of the calculated $`\{u_i^q,v_i^q\}`$ amplitudes onto the condensate. Convergence is reached when the change in $`_i|\stackrel{~}{n}_i|^2`$ from one cycle to the next is smaller than a specified tolerance.

The parameters chosen for the numerical calculations were $`\mathrm{\Omega }=0.0015E_R`$, with $`E_R`$ the one photon recoil energy, which for the case of a Rubidium condensate corresponds to a trap frequency of approximately 90 Hz. We used a total number of 1000 atoms, $`N=1000,`$ and set $`VN=1.0E_R`$. $`J`$ was varied to achieve a range of $`V_{eff}=V/J`$ between 0.01 and 312. The range was chosen based on a local mean field approach , which for our parameters estimates the transition region between $`V_{eff}640`$ (at the center where the local filling factor is approximately 80) and $`V_{eff}12`$ (at the wings).

The results of the numerical calculations are summarized in Figs. 3 to 8. In Fig. 3 we plot the evolution of the density profile (black boxes), the condensate population (triangles) and the on-site depletion (empty diamonds) as $`V_{eff\text{ }}`$ is increased. In the plots we also show, for comparison purposes, the ground state density profile for $`J=0`$ (empty boxes). This has the advantage that can be calculated exactly from the Hamiltonian. In general we observe the reduction of the condensate population and thus the increment of the depletion as the interaction strength is increased. When the system is in the superfluid regime most of the atoms are in the condensate but as $`J`$ is decreased the depletion of the condensate becomes very important.

For the chosen parameters, the density profile has a parabolic shape reflecting the confining potential. By comparing the evolution of the density as $`J`$ is decreased with the exact solution at $`J=0`$, we can crudely estimate the validity of the HFB-Popov calculations. The density evolves from a Gaussian type (see plots for $`V_{eff\text{ }}=0.01`$ and $`0.09`$) with smooth edges towards a Thomas-Fermi profile with sharp edges adjusting its shape to the $`J=0`$ profile. We can appreciate that around $`V_{eff\text{ }}=3`$ both profiles are almost equal. For lower values of $`J`$ the HFB-Popov density starts to differ from the $`J=0`$ one, even though the system is closer to the $`J=0`$ limit. We can say that beyond this point higher order correlations, neglected by the theory, begin to be important. The departure of the HFB-Popov density profile from the $`J=0`$ one as J is decreased begins at the edges (see panel corresponding to $`V_{eff}=11`$ and $`100`$). This is something expected if we look at the on-site depletion. For such values of $`V_{eff}`$ the local depletion at the wings corresponds to a considerable percentage of the condensate populations, and thus the validity of the HFB-Popov assumptions starts to be dubious. The homogeneous results shown in the previous section corroborate our present statements for the confined system. For the smallest filling factor (see Fig.1) the differences between the homogeneous HFB-Popov calculations and the exact solutions become important for values of $`V_{eff\text{ }}`$ greater than 20. For higher values of $`V_{eff}`$, see plot for $`V_{eff\text{ }}=312`$, the HFB-Popov density predictions differs from the $`J=0`$ solution even at the central wells. At this point the failure of the method is clear and a fully quantal method is required.

The HFB-Popov quasiparticle spectrum is shown in Fig.4. It can be observed how the lower energy eigenvalues evolve from a linear non degenerated spectrum to an almost degenerated one as J is decreased. It is worth it to mention that the small energy difference between the ground and first excited states for high values of $`V_{eff\text{ }}`$ makes the numerical solution very unstable in the sense that it was very easy to jump to an excited state when solving for the condensate wave function. The decrement in the energy spacing predicted by the HFB-Popov theory as the system approaches the transition is very useful to keep in mind for the experimental realization of the Mott transition. As the optical lattice depth is ramped up the adiabaticity criteria is harder to fulfill.

In Fig. 5 we plot the results for the number fluctuations found numerically using the inhomogeneous HFB-Popov approach. The number fluctuations profile reflects the condensate profile. We also show the number fluctuations evaluated by using a local density approximation (empty boxes). The latter was calculated by substituting in the number fluctuations expression (Eq. 49) the $`\{u^q,v^q\}`$ amplitudes found for the homogeneous system (Eqs. 51 and 52), but replacing the condensate density in each lattice site by the one found numerically for the trapped system (see Fig. 3). The complete agreement between the two approaches justifies the validity of the local density approximation for the estimations of local quantities in confined systems. Based on this agreement and the results for the homogeneous system shown in the previous section, we expect that the inhomogeneous HFB-Popov results for squeezing also agrees with the exact solution right up to the transition.

### 4.3 Quasi-momentum distribution in the inhomogeneous lattice

The quasi-momentum distribution of the atoms released from the lattice is important because it is one of the most easily accessible quantities to the experiments. The quasi-momentum distribution function $`n_q`$ is defined as

$$n_q=\underset{i,j}{}e^{\mathrm{i}q(ij)a}a_i^{}a_j,$$
(61)
where the quasi-momentum $`q`$ can assume discrete values which are integer multiples of $`\frac{2\pi }{Ia}`$, $`a`$ is the lattice spacing. In Fig. 6 we present the quasi-momentum distribution for the same parameters used in the previous section. The distribution for the two lowest values of $`V_{eff}`$ corresponds to the one that characterizes an uncorrelated superfluid phase with a narrow peak at small quasi-momenta. As the hopping rate is decreased we observe that the sharpness of the central peak decreases and the distribution extends towards large quasi-momenta. It is interesting to note the appearance of a small peak between $`q=0.5`$ and $`1`$ which is most noticeable for the $`V_{eff}=3`$ case. This agrees with the results found in  where they solve numerically the Bose-Hubbard Hamiltonian by using Monte Carlo simulations. We attribute the origin of the small peak to the depletion of the condensate at the wings. For the parameters when the small peak is present, the most important contribution to the quasi-momentum distribution still comes from the condensate atoms. The step function like shape of the condensate profile causes an oscillatory $`|\mathrm{sin}(x)/x|`$ shape of the quasi-momentum distribution. As the lattice depth is increased the hopping becomes energetically costly, the long-range order starts to decrease and the Fourier spectrum becomes broader.

In Fig. 7 we plot the first order on site superfluid fraction $`f_{si}^{(1)}`$ which was defined in Eq.(48). The curves corresponding to $`V_{eff}=0.0111`$, which are in the regime where the HFB-Popov is expected to be valid, depict how as $`V_{eff}`$ is increased the superfluid profile decreases faster at the wings and at the center but no major change is observed in the middle section. The evolution of the on-site superfluidity as the interaction strength is increased, exhibiting a domain localized decrement instead of a global one, is in agreement with the development of uncompressible regions surrounded by superfluid rings predicted for trapped systems  as the transition is approached.

Fig. 8 shows the first order and total superfluid fraction and also the second order superfluid fraction as a function of $`V_{eff}`$. Different from the translationally invariant case, the second order contribution calculated in the HFB-Popov approach doesn’t vanish for the inhomogeneous system. The rapid decrement of the superfluid fraction observed after $`Log(V_{eff})1.2`$ is a signature that the system is entering a highly correlated regime. Beyond this point higher order correlations neglected in the HFB-Popov approach become crucial and a more accurate approach is required.

The Mott transition is a quantum phase transition and as all critical phenomena its behavior depends strongly on the dimensionality of the system. In the present analysis, due to computational limitations, we considered one dimensional systems. Experimentally, the Mott transition has been achieved  in a 3 dimensional lattice with filling factors between 1 and 3. Even though the HFB-Popov approach fails to describe the strong coupling regime for the one dimensional systems we considered in the present paper, we showed how the method is incredibly powerful in describing most of its characteristic features as they are driven from the superfluid regime towards the transition. We expect the HFB-Popov method to give a better description of the transition as the dimensionality of the system is increased and therefore to be a good model in an experimental situation.

As shown in previous studies ,  the Mott transition in a d-dimensional homogeneous system has two different critical behaviors: one (d+1) XY- like, for systems with fixed integer density as the interaction strength is changed, and one mean field-like exhibited when the transition is induced by changing the density. Different from the homogeneous case where the Mott transition is characterized by the global offset of the superfluidity, for confined systems, commensuration is only well defined locally. The inhomogeneity introduced by the confined potential allows the existence of extended Mott domains (above a critical interaction strength) surrounded by superfluid ones , thus the total superfluid fraction doesn’t vanishes in the Mott regime. This issue, together with the fact that the finite length scale introduced by the trap suppresses the long wave fluctuations which are responsible for destroying the mean field <sup>1</sup><sup>1</sup>1One obvious consequence of this is that BEC is possible in one and two dimensions in a trap whereas in the homogeneous, thermodynamic limit it can not occur in fewer than three dimensions, make us believe the critical behavior in confined systems to be more mean-field like. Because the critical dimension for the latter type of transition is two , , we expect that for trapped systems in $`d=3`$, the range of validity of the HFB-Popov extends closer to the transition.

## 5 Summary

We have developed in this article a Bogoliubov method for describing the approach of a condensate loaded in an optical lattice towards the Mott transition. We have shown that this method can be used to predict the relevant physical quantities over a useful range. We have also shown how it gives a powerful insight into the way quantum depletion reduces the long range order and the superfluid fraction.

## Acknowledgments

This work was supported in part by US National Science Foundation grants PHY-0100634 and PHY-0100767, the United Kingdom’s Engineering and Physical Sciences Research Council, the Cold Quantum Gases Network Research Training Network, and the Advanced Research and Development activity.


format: text

text:
# Tachyon Dynamics and Brane Cosmology

## Acknowledgments

It is my pleasure to thank Ira Wasserman and Henry Tye for enjoyable collaborations. This work was supported in part by the Department of Energy.


text:
# Percolating cluster of center vortices and confinement

## 1 INTRODUCTION

Center vortices are believed to play an important role in confining gauge theories: they disorder the gauge configurations and hence make the Wilson loop decay with the area law. This idea is quite old , but its investigation by lattice methods began recently .

We study a very simple, though non trivial, gauge theory, namely $`_2`$ gauge theory in three space–time dimensions: it turns out to be an ideal laboratory to to test the relevance of center vortices to the confinement mechanism; in fact, in this theory, the gauge group coincides with its center, so the identification of center vortices can be performed without resorting to any gauge-fixing procedure. Moreover, the rather small size of configuration space enables us to get high precision numerical results from Monte Carlo simulations.

We focus our attention on the effect of the size of clusters of center vortices on confinement: the hypothesis we want to test, first proposed in Ref. , is that confinement in this model is due to the existence of an infinite cluster of center vortices, i.e. a connected component, in the graph defined in the dual lattice by all the center vortices, whose size scales linearly with the lattice volume. Only such a giant component in the graph can disorder the gauge configurations enough to produce the area-law decay of the Wilson loop.

In order to shed light on the role of large clusters of center vortices, we first modify each configuration of a Monte Carlo-generated ensemble in the confined phase by removing the largest cluster of center vortices, and show that the ensemble thus obtained does not confine. On the other hand, if we remove all of the small clusters of center vortices, and leave the largest one only, confinement is preserved, albeit with a string tension significantly smaller than the original one. Remarkably, also the string corrections due to the quantum fluctuations of the confining flux tube are preserved by this transformation.

## 2 THE MODEL

The $`_2`$ gauge model is defined by the action

$$S(\beta )=\beta \underset{\mathrm{}}{}\sigma _{\begin{array}{c}\\ \mathrm{}\end{array}},\sigma _{\begin{array}{c}\\ \mathrm{}\end{array}}=\underset{\mathrm{}\mathrm{}}{}\sigma _{\mathrm{}}$$
(1)
where the sum is extended to all plaquettes of a cubic lattice; the fundamental degrees of freedom are $`_2`$ variables $`\sigma `$ defined on links $`\mathrm{}`$, and each plaquette contributes the product of its links to the action. The model enjoys a $`_2`$ gauge invariance: any local transformation consisting in flipping all of the $`\sigma `$’s defined on the links that meet at a given site, leaves the action invariant.

We can construct center vortices by assigning a vortex line in the dual lattice to each frustrated plaquette in the direct lattice. The resulting graph of center vortices in the dual lattice is gauge–invariant, and has an even coordination number.

Such graph is in general made of many connected components: the clusters of center vortices. The value of a Wilson loop in a given configuration is $`\pm 1`$ according to the number, modulo 2, of vortex lines that are linked to the loop. A plausibility argument for the crucial role played by the infinite cluster of center vortices is the following (see also Ref. ): the only clusters contributing to the Wilson loop $`W(C)`$ are those linked to $`C`$; so, in the limit of large enough Wilson loops, clusters of finite size $`R`$ can contribute to $`W(C)`$ only when they are near $`C`$, and their number grows linearly with the length of the loop: therefore, $`W(C)`$ decays with a perimeter law, and the theory is deconfined. On the other hand, if there exists an infinite cluster of center vortices, then the number of vortex lines linked to the loop grows with its area, and the system is in the confining regime. It is easy to verify numerically that in the confining phase of $`_2`$ gauge system there is no ambiguity in finding a cluster of center vortices whose size scales linearly with the lattice volume for large enough lattices; conversely, in the deconfined phase the density of the largest cluster decreases rapidly with the volume.

## 3 ROLE OF THE LARGE CLUSTER OF CENTER VORTICES

To study the relation between the infinite cluster and the string tension value, we chose to simulate the model at $`\beta =0.74883`$, which is well inside the scaling region, and for which the value of the string tension is known with high precision : $`\sigma =0.01473(10)`$.

First, we verified that, in this confining regime, in the dual lattice configurations there exists a large cluster of center vortices, whose size grows linearly with the lattice volume.

Since VEV’s of rectangular Wilson loops of sizes $`R\times T`$ in the confining regime of a gauge theory are expected to behave as:

$$W(R,T)=e^{\sigma RT+p(R+T)+k}\left[\frac{\eta \left(i\frac{T}{R}\right)}{\sqrt{R}}\right]^{\frac{d2}{2}}$$
(2)
(where the term in square brackets is associated with quantum fluctuations of the colour flux tube), an efficient method  to extract the string tension $`\sigma `$ from Wilson loop data generated by Monte Carlo simulations consists in considering loops of the same perimeter, but with different areas; we can define the following quantity:

$$\sigma _{\mathrm{eff}}(L,n)\frac{1}{n^2}\mathrm{log}\left[\frac{W(L+n,Ln)}{W(L,L)F(n/L)}\right]$$
(3)
where:

$$F(t)\left[\frac{\eta (i)\sqrt{1t}}{\eta \left(i\frac{1+t}{1t}\right)}\right]^{1/2}$$
(4)
The estimator $`\sigma _{\mathrm{eff}}`$ approaches the string tension for large $`L`$, when it is evaluated by data obtained from the Monte Carlo ensemble of configurations.

Now, if we erase the largest cluster of center vortices (in the dual lattice) in each of these configurations, the string tension estimator $`\sigma _{\mathrm{eff}}`$ approaches zero for large $`L`$: this means that removing the “infinite” cluster leaves us with non-confining configurations.

On the other hand, if we erase all clusters but the largest one, we still obtain confining configurations, but the string tension evaluated on these new configurations turns out to be significantly smaller than its value in the original Monte Carlo configurations.

Another very interesting result (see  for details) is that the corrections due to the flux tube fluctuations survive the elimination of the small clusters.

These results may be compared with those obtained in other gauge theories .

## 4 CONCLUSIONS

Our numerical results show that:

* The largest center vortex in the dual lattice is responsible for confinement, but the string tension measured from configurations in which all other clusters have been removed does not reproduce the full string tension of the original theory: small clusters of vortices are unable by themselves to disorder the system enough to produce confinement, but they do give a finite contribution to the string tension of the full theory.
* Quantum fluctuations of the flux tube survive the elimination of the small clusters: after deletion of all the small clusters, the Wilson loop VEV shows the same shape dependence as in the full theory, which can be explained as originating by the fluctuations of a free bosonic string.


text:
# Effective sigma models and lattice Ward identities

## 1 Introduction

Confinement in pure Yang–Mills theory is still a theoretical challenge. The problem actually has two faces. On the one hand, there is confinement of static external sources in the fundamental representation which manifests itself through the appearance of a linear potential (nonzero string tension). On the other hand, there should be gluon confinement implying a finite range of the gluonic interactions, i.e. a mass gap. How the two faces are related is largely unclear at the moment.

Recently, Faddeev and Niemi (FN) have suggested that the infrared dynamics of glue might be described by the following low–energy effective action ,

$$S_{\mathrm{FN}}=d^4x\left[m^2(_\mu 𝒏)^2+\frac{1}{e^2}H_{\mu \nu }H^{\mu \nu }\right].$$
(1)
Here, $`𝒏`$ is a unit vector field with values on $`S^2`$, $`𝒏^2n^an^a=1`$, $`a=1,2,3`$; $`m`$ is a dimensionful and $`e`$ a dimensionless coupling constant. The FN ‘field strength’ is defined as

$$H_{\mu \nu }𝒏_\mu 𝒏\times _\nu 𝒏.$$
(2)
Faddeev and Niemi argued that (1) “is the unique local and Lorentz–invariant action for the unit vector $`𝒏`$ which is at most quadratic in time derivatives so that it admits a Hamiltonian interpretation and involves all such terms that are either relevant or marginal in the infrared limit”.

It has been shown that $`S_{\mathrm{FN}}`$ supports string–like knot solitons , characterized by a topological charge which equals the Hopf index of the map $`𝒏:S^3S^2`$. Here, $`𝒏`$ is supposed to be static and approaches a uniform limit at spatial infinity, $`𝒏_{\mathrm{}}=𝒆_z`$. In analogy with the Skyrme model, the $`H^2`$ term is needed for stabilization. The knot solitons can possibly be identified with closed gluonic flux tubes and are thus conjectured to correspond to glueballs. For a rewriting in terms of curvature–free $`SU(2)`$ gauge fields and the corresponding reinterpretation of $`S_{\mathrm{FN}}`$ we refer to .

In order for the model to really make sense, however, the following problems have to be solved. First of all, neither the interpretation of $`𝒏`$ nor its relation to Yang–Mills theory have been fully clarified. An analytic derivation of the FN action requires

* an appropriate change of variables, $`A(𝒏,X)`$, relating the Yang–Mills potential $`A`$ to $`𝒏`$ and some remainder $`X`$
* the functional integration over $`X`$ to arrive at an effective action $`S_{\mathrm{eff}}`$ for the $`𝒏`$-field.

Some progress in this direction has been made  on the basis of the Manton–Cho decomposition ,

$$𝑨_\mu =C_\mu 𝒏𝒏\times _\mu 𝒏+𝑾_\mu ,$$
(3)
where $`C`$ is an Abelian connection and $`𝒏𝑾_\mu =0`$. Nevertheless, it is fair to say that there are no conclusive results up to now.

Second, there is no reason why in a low–energy effective action for the $`𝒏`$–fields both operators in the FN ‘Skyrme term’, which can be rewritten as

$$H^2=(_\mu 𝒏_\mu 𝒏)^2(_\mu 𝒏_\nu 𝒏)^2,$$
(4)
should have the same coupling. Third, and conceptually most important, $`S_{\mathrm{FN}}`$ has the same spontaneous symmetry breaking pattern as the nonlinear $`\sigma `$–model, $`SU(2)U(1)`$. Hence, it should admit two Goldstone bosons and one expects to find no mass gap. In order to exclude these unwanted massless modes we have recently suggested to break the global $`SU(2)`$ explicitly , an idea that has subsequently also been adopted by Faddeev and Niemi .

In what follows the FN hypothesis will be tested on the lattice. To avoid the appearance of Goldstone bosons we allow for explicit symmetry–breaking terms.

## 2 Generating an $`𝑺𝑼\mathbf{(}\mathrm{𝟐}\mathbf{)}`$ lattice ensemble of $`𝒏`$–fields

The conceptual problem to be solved in the first place is to obtain a reasonable ensemble of $`𝒏`$–fields. The (lattice version of the) decomposition (3) is of no help: it assumes some particular choice of $`𝒏`$ on which the decomposition is then based. One way of defining an $`𝒏`$–field is via Abelian gauge fixing, originally introduced by ‘t Hooft . A prominent example in this class of gauges is the maximally Abelian gauge (MAG) which is obtained via maximizing the functional

$$F_{\mathrm{MAG}}[U;g]\underset{x,\mu }{}\text{tr}\left(\tau _3^gU_{x,\mu }\tau _3^gU_{x,\mu }^{}\right)\underset{x,\mu }{}\text{tr}\left(n_xU_{x,\mu }n_{x+\mu }U_{x,\mu }^{}\right)\stackrel{~}{F}_{\mathrm{MAG}}[U;n],$$
(5)
with respect to the gauge transformation $`g`$. The maximizing $`g`$ then defines the $`𝒏`$–field according to

$`n_xg_x^{}\tau _3g_x𝒏_x𝝉.`$ (6)
Instead of maximizing $`F_{\mathrm{MAG}}`$ with respect to $`g`$ one can equivalently maximize $`\stackrel{~}{F}_{\mathrm{MAG}}`$ with respect to $`n`$  which results in the condition

$$\mathrm{}[U]n_x\lambda _xn_x.$$
(7)
Here, $`\mathrm{}[U]`$ denotes the covariant Laplacian in the adjoint representation (see App. A), while $`\lambda _x`$ is a Lagrange multiplier imposing that $`𝒏_x`$ is normalized to unity, locally at each lattice site $`x`$. In principle, (7) can be solved for the field $`𝒏`$ associated with the background $`U`$. However, as this background is distributed randomly along its orbit we will in turn obtain a random ensemble of $`𝒏`$–fields characterized by the two–point function

$$G_{xy}^{ab}n_x^an_y^b=\frac{1}{3}\delta ^{ab}\delta _{xy}.$$
(8)
Thus, nontrivial correlations are absent. Fig. 1 shows that this is indeed what one gets in a typical Monte Carlo run.

One way out of this problem is to follow the continuum approach of  which starts out with a covariant gauge fixing. After having generated $`SU(2)`$ lattice configurations using the standard Wilson action we therefore fix to lattice Landau gauge (LLG). The latter is defined by maximizing the functional

$$F_{\mathrm{LLG}}\underset{x,\mu }{}\text{tr}^\mathrm{\Omega }U_{x,\mu },$$
(9)
with respect to the gauge transformation $`\mathrm{\Omega }`$. In this way we impose some ‘preconditioning’  which (i) eliminates the randomness in our Yang–Mills ensemble and (ii) leaves a residual global $`SU(2)`$–symmetry. The Landau gauge configurations are then plugged into the MAG functional (5) which subsequently is maximized with respect to $`g`$. The gauge transformation $`g`$ obtained this way determines $`𝒏`$ according to (6). One may say that $`g`$ (and hence $`𝒏`$) measure the gauge–invariant (!) distance between the LLG and MAG gauge slices (see Fig. 2).

In App. B we show that LLG and MAG are ‘close’ to each other. Therefore, the maximizing $`g`$ is on average close to unity, hence, on average, $`𝒏`$ will be aligned in the positive 3–direction. In this way we have explicitly broken the global $`SU(2)`$ down to a global $`U(1)`$.

All computations have been done on a $`16^4`$–lattice with Wilson coupling $`\beta =2.35`$, lattice spacing $`0.13\mathrm{fm}`$ and periodic boundary conditions. For the LLG we used Fourier accelerated steepest descent  (see Fig. 3). The MAG was achieved using two independent algorithms, one (AI) using iterations based on elementary geometric manipulations (including overrelaxation steps), the other (AII) being analogous to LLG fixing (see Fig. 4).

## 3 Numerical results

As expected, we observe a non–vanishing expectation value of the field in the 3–direction, a ’magnetization’ $`𝔐`$ defined through $`n^a=𝔐\delta ^{a3}`$. Thus, the global symmetry is indeed broken explicitely according to the pattern $`SU(2)U(1)`$. We demonstrate this by exhibiting the angular distribution of the $`𝒏`$–field on its target space $`S^2`$ in Fig. 5. The azimuthal angle $`\varphi `$ is equally distributed, while the distribution of the polar angle $`\theta `$ has a maximum near $`\pi /2`$ corresponding to the north pole, $`𝒏=(0,0,1)`$.

Explicit symmetry breaking also shows up in the behavior of the two–point functions (Fig. 6). The longitudinal correlator, $`G_x^{}n_x^3n_0^3n^3n^3=𝔐^2`$, exhibits clustering for large distances, the plateau being given by the magnetization (squared).

The transverse correlation function (of the would–be Goldstone bosons)

$`G_x^{}G_{x0}^{}{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{2}{}}}n_x^in_0^i,`$ (10)
decays exponentially as shown in Fig. 7. This means that there is a nonvanishing mass gap $`M`$ whose value can be obtained by a fit to a $`\mathrm{cosh}`$–function (see Fig. 7).

The numerical values of the observables, $`𝔐`$, $`M`$ and the transverse susceptibility,

$$\chi ^{}\underset{x}{}G_x^{},$$
(11)
which all can be derived from the two–point functions, are summarized in Table 1 for both algorithms.

The disagreement between AI and AII is statistically significant. We attribute it to the ubiquitous Gribov problem  (for Abelian gauges, see ). On the lattice, this is the statement that maximizing gauge fixing functions like $`F_{\mathrm{MAG}}`$ or $`F_{\mathrm{LLG}}`$ is equivalent to a spin–glass problem with an enormous number of degenerate extrema. This implies that the algorithms AI and AII will almost certainly end up in different local maxima, which explains the difference between rows one and two in Table 1.

As shown in the last column of Table 1, the numerical results for the mass gap $`M`$ lead to a value of about 1 GeV in physical units.

$`n^3`$ is a local functional of the $`n^i`$, $`n^3=(1n^in^i)^{1/2}`$. Thus, one expects the same exponential decay for the longitudinal correlator $`G^{}`$. This can be confirmed with a numerical value for the mass gap of $`M=0.66a^1`$.

To improve statistics, we have calculated the time–slice correlator,

$$C^{}(t)L^3\underset{𝒙}{}G_{𝒙,t}^{}.$$
(12)
In the continuum, for purely exponential decay of $`G^{}`$, this would become proportional to a modified Bessel function $`K_2`$. An associated fit works very well as is shown in Fig. 8.

Fitting the time–slice correlator according to Fig. 8, we obtain for the mass gap

$$aM=0.642\text{i.e.}M=0.97\text{GeV}.$$
(13)
This is the value with the smallest statistical errors.

The mass gap obtained differs significantly from the $`SU(2)`$ mass gap, $`M_{SU(2)}`$ 1.5 GeV, obtained directly from a Wilson ensemble with $`\beta =2.4`$ . We believe that the difference is due to the highly nonlocal relation between the original Yang–Mills degrees of freedom (the link variables) and the color spin $`𝒏`$. After all, we have implicitly solved the partial differential equation (7) with link variables $`U`$ in Landau gauge entering the adjoint Laplacian. The solution $`𝒏`$ will clearly be a nonlocal functional of these $`U`$’s. Consequently, we cannot expect that the exponential decay of $`G^{}[n]`$ will be governed by the lowest excitation of the $`U`$–ensemble.

## 4 Effective action and Schwinger–Dyson equations

At this point it is natural to ask whether there is an effective action $`S_{\mathrm{eff}}[n]`$ that reproduces the distribution of $`𝒏`$–fields leading to the results of the previous section.

At low energies, it should make sense to employ an ansatz in terms of a derivative expansion,

$$S_{\mathrm{eff}}=\underset{j}{}\lambda _jS_j[𝒏]+\underset{j}{}\lambda _j^{}S_j^{}[𝒏,𝒉],$$
(14)
with $`O(3)`$ invariant operators $`S_j`$ and noninvariant operators $`S_j^{}`$, which are ordered by increasing mass dimension. Up to dimension four, one has the symmetric terms,

$$\begin{array}{ccc}S_1=(𝒏,\mathrm{}𝒏),\hfill & & S_2=(𝒏,\mathrm{}^2𝒏),\hfill \\ S_3=(𝒏\mathrm{}𝒏,𝒏\mathrm{}𝒏),\hfill & & S_4=(𝒏_\mu ^{}_\nu 𝒏,𝒏_\mu ^{}_\nu 𝒏),\hfill \end{array}$$
(15)
and the symmetry–breaking terms including a unit vector ‘source field’ $`𝒉`$  (which can be thought of as the direction of an external magnetic field),

$`S_1^{}=(𝒏,𝒉),S_2^{}=(𝒏𝒉,𝒏𝒉),S_3^{}=(𝒏\mathrm{}𝒏,𝒏𝒉).`$ (16)
In the above, we have introduced the scalar products

$$(f,g)\underset{x}{}f_xg_x,𝒖𝒗u^av^a,$$
(17)
and the usual lattice Laplacian $`\mathrm{}`$ (see App. A).

Note that the $`𝒏`$–field configurations are classified by the Hopf invariant irrespective of the particular form of the (effective) action. This, together with the usual scaling arguments, shows that the action (14) with the operators (15) and (16) should still support classical knot soliton solutions. Our ansatz thus does not exclude this important feature.

The couplings in (14) can be determined by inverse Monte Carlo techniques. The notion is suggestive: instead of creating an ensemble from a given action, one wants to compute a (truncated) action which gives rise to the given ensemble of $`𝒏`$–fields. A particular approach uses the Schwinger–Dyson equations . These represent an overdetermined linear system which can be used to solve for the couplings in terms of correlation functions. The latter are nothing but the coefficients of the linear system.

For an unconstrained scalar field $`\varphi `$, the Schwinger–Dyson equations follow from translational invariance of the functional measure, implying

$$0=𝒟\varphi _x\left\{F[\varphi ]\mathrm{exp}(S[\varphi ])\right\},$$
(18)
where $`_xi\delta /\delta \varphi _x`$ is the (functional) momentum operator, and $`F`$ an arbitrary functional111Usually one chooses $`F[\varphi ]=\varphi (x_1)\mathrm{}\varphi (x_k)`$. of the field $`\varphi `$. For a constrained field like $`𝒏`$ with a curved target space things are slightly more subtle . There is, however, a rather elegant way to derive the Schwinger–Dyson equations if one exploits the isometries of the target space $`S^2`$ . The target space measure,

$$𝒟𝒏=\underset{x}{}d𝒏_x\delta \left(𝒏_x^21\right),$$
(19)
is obviously rotationally invariant, i.e. under $`𝒏R𝒏`$, $`RO(3)`$. This implies the modified Schwinger–Dyson identity

$$𝒟𝒏𝑳_x\left\{F[𝒏]\mathrm{exp}\left(S_{\mathrm{eff}}[𝒏,𝒉]\right)\right\}=0,$$
(20)
where $`𝑳_x`$ denotes the angular momentum operator (at lattice site $`x`$),

$$i𝑳_x=𝒏_x\times \frac{}{𝒏_x}\text{or}iL_x^aϵ^{abc}n_x^b\frac{\delta }{\delta n_x^c}.$$
(21)
In shorthand–notation, (20) can be rewritten as

$$𝑳_xF[𝒏]F[𝒏]𝑳_xS_{\mathrm{eff}}[𝒏]=0.$$
(22)
These exact identities can be used to determine the unknown couplings $`\lambda _j`$. To this end one chooses a set of field monomials $`F_i[𝒏]`$ and plugs them into (22) together with the form (14) of the action. This yields the local linear system

$$\underset{j}{}F_i𝑳_xS_j\lambda _j+\underset{k}{}F_i𝑳_xS_k^{}\lambda _k^{}=𝑳_xF_i,$$
(23)
which, in principle, can be solved numerically, for instance by least–square methods. The identities obtained so far hold for arbitrary actions $`S_{\mathrm{eff}}[𝒏]`$. In particular, we have not made use of any symmetries. Taking the latter into account will lead to Ward identities.

Let us specialize to our lattice effective action (14). It is a sum of a symmetric part $`S`$ containing the terms (15) and an asymmetric part $`S^{}`$ containing the terms (16),

$$S_{\mathrm{eff}}=S[𝒏]+S^{}[𝒏,𝒉].$$
(24)
Due to the invariance of $`S`$ under global $`O(3)`$ rotations it is an $`O(3)`$–singlet and hence annihilated by the total angular momentum,

$$𝑳S=0,𝑳=\underset{x}{}𝑳_x,$$
(25)
such that $`𝑳S_{\mathrm{eff}}=𝑳S^{}`$. Thus, summing over all lattice sites $`x`$ in (22) yields the (broken) Ward identity,

$$𝑳F[𝒏]F[𝒏]𝑳S^{}[𝒏,𝒉]=0,$$
(26)
where the second terms contains the infinitesimal change of the non-invariant part $`S^{}`$ of the effective action under rotations of $`𝒏`$. Note that the coupling constants $`\lambda _j`$ of the $`O(3)`$–symmetric operators $`S_j`$ have disappeared in the Ward identity (26) so that only the symmetry–breaking couplings $`\lambda ^{}`$ are present. We have collected the explicit lattice Schwinger–Dyson and Ward identities used in our simulations in App. C. As the former are local relations, they naturally contain more information than the global Ward identities. In particular, one does have access to length scales.

## 5 Comparing Yang–Mills and FN ensembles

### 5.1 Leading–order ansatz

To leading order (LO) in the derivative expansion we have a standard nonlinear sigma model with symmetry–breaking term,

$$S_{\mathrm{eff}}=\underset{x}{}(\lambda 𝒏_x\mathrm{}𝒏_x+\lambda ^{}𝒏_x𝒉),𝒉𝒆_z.$$
(27)
Inverse Monte Carlo amounts to determining the couplings $`\lambda `$ and $`\lambda ^{}`$ such that the probability distribution associated with the LO action (27) fits the observables of the Yang–Mills ensemble of $`𝒏`$–fields222Throughout this section, we refer to algorithm AI.. The associated Schwinger–Dyson equation (23), with $`F[𝒏]=n_x^a`$, can be written as

$$\lambda H_{xy}+\lambda ^{}G_{xy}^{}=𝔐\delta _{xy},$$
(28)
where $`H`$ denotes the (antisymmetrized) two–point function of $`n^i`$ and $`n^i\mathrm{}n^3`$,

$$H_{xy}n_x^in_y^i\mathrm{}n_y^3n_x^in_y^3\mathrm{}n_y^in_x^in_y^{[i}\mathrm{}n_y^{3]}.$$
(29)
To analyse (28) we define a ‘reduced’ two–point function $`h_{xy}`$ and magnetization $`\mu `$,

$$h_{xy}H_{xy}/G_{xy}^{},\mu 𝔐/G_{xx}^{},$$
(30)
and rewrite (28) as the inhomogeneous system (using translational invariance to replace $`xyx`$),

$`\lambda h_x+\lambda ^{}`$ $`=`$ $`0,x=1,\mathrm{},8,`$ (31)
$`\lambda h_0+\lambda ^{}`$ $`=`$ $`\mu .`$ (32)
The solution is found to be

$`\lambda `$ $`=`$ $`{\displaystyle \frac{1}{h_xh_0}}\mu ,`$ (33)
$`\lambda ^{}`$ $`=`$ $`{\displaystyle \frac{h_x}{h_xh_0}}\mu ,`$ (34)
with the numerical boundary value given by $`h_0=0.1410`$ (cf. Fig. 9). Clearly, the system (31), (32) is overdetermined (nine equations for two unknowns). This is reflected in the fact that $`\lambda `$ and $`\lambda ^{}`$ in (33) and (34) depend on the lattice distance $`x`$ via $`h_x`$. If the Yang–Mills ensemble were exactly described by the LO action (27), there would be no such $`x`$–dependence. Rather, for any $`x=1,\mathrm{},8`$, we would have the same values for $`\lambda `$ and $`\lambda ^{}`$, respectively. Thus, to test the quality of the LO ansatz, we divide (34) by (33) showing explicitly that $`h_x`$ should be constant,

$$h_x=\frac{\lambda ^{}}{\lambda }\kappa ^{}=const.,x0.$$
(35)
Fig. 9 shows that this is not the case.

Therefore, a minimal sigma model with symmetry–breaking term does not yield a good representation of our Yang–Mills ensemble of $`𝒏`$–fields. If we nevertheless insist on the LO description, we have to ‘fit’ $`h_x`$ by a horizontal line so that the numerical determination of the couplings via (33) and (34) is beset by large errors,

$`\lambda `$ $`=`$ $`1.41\pm 5.25,`$ (36)
$`\lambda ^{}`$ $`=`$ $`1.33\pm 0.74.`$ (37)
Obviously, $`\lambda `$ (including its sign) remains essentially undetermined. For $`\lambda ^{}`$ the situation is slightly better.

In order to assess the errors it is worthwhile to check whether our numerical accuracy is sufficient to really validate the Schwinger–Dyson identity (28) for the LO action (27) on the lattice. To this end we have simulated (27) with a combination of Metropolis and cluster algorithms producing a number of 150 configurations using the central values (36) and (37) as the input couplings. The result for $`h_x`$ in the LO ensemble is presented in Fig. 10.

It is reassuring to note that the simulation of the minimal sigma model reproduces the input value $`\kappa ^{}=0.943`$ very well (for $`x0`$), the error being of the order of one percent. The prediction (35) thus can be verified with high accuracy for the LO action (27). We conclude that inverse Monte Carlo works quite well when applied to the minimal $`\sigma `$–model.

The discrepancy between the LO and Yang–Mills ensembles can be further visualized by looking at the susceptibility. For the action (27) and the choice $`F[𝒏]=n_x^a`$, the Ward identity (26) assumes the simple form

$$\chi ^{}=𝔐/\lambda ^{}.$$
(38)
A consistency check is provided by noting that this can directly be obtained by summing (28) over $`x`$. Plugging in the magnetization from Table 1 and $`\lambda ^{}`$ from (37) we find

$$\chi ^{}=0.33\pm 0.18,$$
(39)
This is way off the Yang–Mills value of 92.6 displayed in Table 1. For magnetization and mass gap the simulation of the LO ensemble yields the values

$$𝔐=0.93,M=1.5,$$
(40)
which are both larger than the Yang–Mills values of Table 1.

The discussion of this subsection thus shows quite clearly that more operators will have to be included in order to possibly make inverse Monte Carlo work reasonably well.

### 5.2 FN action with symmetry–breaking term

In this subsection we consider the FN action (1) with a LO symmetry–breaking term,

$$S_{\mathrm{eff}}=\underset{x}{}\left\{\lambda 𝒏_x\mathrm{}𝒏_x+\lambda _{\text{FN}}\left[(𝒏\mathrm{}𝒏)^2(𝒏_\mu ^{}_\nu 𝒏)^2\right]+\lambda ^{}𝒏_x𝒉\right\}.$$
(41)
This ansatz does not include all terms of next–to–leading order (NLO) in the derivative expansion. It should be viewed as a minimal modification of the original FN action by adding an explicit symmetry–breaking term to obtain a mass gap.

The Schwinger–Dyson equation generalizing (28) becomes

$$\lambda H_{xy}+\lambda _{\text{FN}}H_{xy}^{\text{FN}}+\lambda ^{}G_{xy}^{}=𝔐\delta _{xy}.$$
(42)
The new two–point function $`H^{\text{FN}}`$ is given by (90). The local identities (42) are to be solved for the three unknown couplings $`\lambda `$, $`\lambda ^{}`$ and $`\lambda _{\text{FN}}`$.

Introducing another reduced two–point function,

$$h_{xy}^{\text{FN}}H_{xy}^{\text{FN}}/G_{xy}^{},$$
(43)
which is plotted in Fig. 11, we obtain, instead of (31) and (32), the (overdetermined) system,

$`\lambda h_x+\lambda _{\text{FN}}h_x^{\text{FN}}+\lambda ^{}`$ $`=`$ $`0,x=1,\mathrm{},8;`$ (44)
$`\lambda h_y+\lambda _{\text{FN}}h_y^{\text{FN}}+\lambda ^{}`$ $`=`$ $`0,y>x;`$ (45)
$`\lambda h_0+\lambda _{\text{FN}}h_0^{\text{FN}}+\lambda ^{}`$ $`=`$ $`\mu ,`$ (46)
to be solved for each pair of lattice distances $`(x,y)`$, $`y>x`$. The number of independent pairs is $`7(7+1)/2=28`$ for lattice extension $`L=16`$. The solutions, labelled by $`x`$ and $`y`$, are

$`\lambda `$ $`=`$ $`{\displaystyle \frac{h_x^{\text{FN}}h_y^{\text{FN}}}{d_{xy}}}\mu ,`$ (47)
$`\lambda ^{}`$ $`=`$ $`{\displaystyle \frac{h_xh_y^{\text{FN}}h_yh_x^{\text{FN}}}{d_{xy}}}\mu ,`$ (48)
$`\lambda _{\text{FN}}`$ $`=`$ $`{\displaystyle \frac{h_xh_y}{d_{xy}}}\mu ,`$ (49)
where we have defined the determinant

$$d_{xy}h_0(h_x^{\text{FN}}h_y^{\text{FN}})h_0^{\text{FN}}(h_xh_y)+h_xh_y^{\text{FN}}h_yh_x^{\text{FN}}.$$
(50)
We have checked that the numerical values for $`d_{xy}`$ are not close to zero so that there is no problem with small denominators in the solutions.

For each pair of lattice distances $`(x,y)`$ we thus have a certain value for any of the three couplings (4749). For each particular coupling those would all agree (within statistical errors) if the NLO action (41) would exactly describe the Yang–Mills ensemble. Again, however, analogous to the LO case, the couplings do vary with lattice distances $`x`$ and $`y`$ as shown in Fig. 12. Numerically, one finds,

$`\lambda `$ $`=`$ $`0.232\pm 0.035,`$ (51)
$`\lambda ^{}`$ $`=`$ $`0.257\pm 0.014,`$ (52)
$`\lambda _{\text{FN}}`$ $`=`$ $`0.0402\pm 0.0004.`$ (53)
Several remarks are in order. First of all, the relative errors, given by the standard deviation from the mean (see Fig. 12), are small compared to the LO ansatz. In particular, the signs of all couplings are fixed. Interestingly, the addition of the FN coupling $`\lambda _{\text{FN}}`$, although small numerically, has a large effect: it reverts the sign of $`\lambda ^{}`$ as compared to (37), implying a negative magnetization. This follows, for instance, from the Ward identity (38), which still holds for the action (41), and the positivity of the susceptibility, hence

$$𝔐=\lambda ^{}\chi ^{}<0,$$
(54)
in contradistinction with the positive Yang–Mills value of Table 1.

To further analyse the result for the couplings, we divide (44) by $`\lambda `$, leading to a linear relation between $`h`$ and $`h^{\text{FN}}`$ (for $`x0`$),

$$h_x=\kappa ^{}\kappa _{\text{FN}}h_x^{\text{FN}},\kappa _{\text{FN}}\lambda _{\text{FN}}/\lambda .$$
(55)
Thus, plotting $`h_x`$ against $`h_x^{\text{FN}}`$ should yield a straight line with intercept $`\kappa ^{}`$ and slope $`\kappa _{\text{FN}}`$. The numerical values (5153) yield

$$\kappa ^{}=1.108\pm 0.228,\kappa _{\text{FN}}=0.173\pm 0.024.$$
(56)
In analogy with the LO case, we have numerically checked the prediction (55) for the NLO action (41) by a Monte Carlo simulation with 150 configurations using the input couplings (5153). Fig. 13 clearly demonstrates the expected linear behavior. A corresponding fit results in

$$\kappa ^{}=1.120,\kappa _{\text{FN}}=0.171,$$
(57)
being consistent with the central values of (56) to within one percent.

We thus conclude that inverse Monte Carlo also works quite well for the NLO ensemble. For the sake of explicit comparison with Fig.s 9 and 11 we display the reduced two–point functions obtained by simulating the NLO action in Fig.s 14 and 15.

As expected, the NLO simulation yields a negative magnetization,

$$𝔐=0.49,$$
(58)
while the mass gap becomes $`M=1.2`$, i.e. slightly larger than the value listed in Table 1.

If the Yang–Mills ensemble has anything to do with the FN one, then plotting $`h`$ vs. $`h^{\text{FN}}`$ (as obtained from Yang–Mills) should also show straight–line behavior, at least approximately. Fig. 16 displays a linear fit to the the Yang–Mills data with parameters

$$\kappa ^{}=1.232,\kappa _{\text{FN}}=0.170.$$
(59)
Again, within error bars, these values are consistent with the preceding analysis (56). Note that reverting the sign of $`\lambda ^{}`$ amounts to reverting the sign of the intercept in Fig. 16. The data points clearly do not support anything like that. On the contrary, it seems that, for small $`h^{\text{FN}}`$ (corresponding to small distances $`x`$, see Fig. 11), the data points deviate from a straight line. Playing around with different fits indicates that $`h`$ rises with a negative power of $`h^{\text{FN}}`$ for small $`h^{\text{FN}}`$ so that in reality there may be no intercept at all.

The behavior of $`\kappa ^{}`$ as a function of lattice distances $`x`$ and $`y`$ may also be investigated. Dividing (48) by (47) we find

$$\kappa ^{}=\lambda ^{}/\lambda =\frac{h_xh_y^{\text{FN}}h_yh_x^{\text{FN}}}{h_x^{\text{FN}}h_y^{\text{FN}}}=1.38\pm 0.29,$$
(60)
where, in the analytic expression, the determinant $`d_{xy}`$ has dropped out. Fig. 17 shows the variation of $`\kappa ^{}`$ with $`x`$ and $`y`$. Again, a different sign for $`\kappa ^{}`$ is completely out of reach.

Following the logic of gradient expansions, one may argue that the effective action (41) is supposed to represent the Yang–Mills ensemble only for large distances. Fig. 16, for instance, seems to indicate that the straight–line fit works particularly well for the last three points to the right which correspond to $`x=6,7,8`$, respectively. In physical units, this amounts to distances $`R`$ larger than six lattice units, i.e. $`R`$ $`\stackrel{>}{}`$ 0.8 fm. Restricting to the analogous data points in Fig. 12, we obtain for the couplings in (41),

$`\lambda `$ $`=`$ $`0.2775\pm 0.026,`$ (61)
$`\lambda ^{}`$ $`=`$ $`0.2661\pm 0.012,`$ (62)
$`\lambda ^{\text{FN}}`$ $`=`$ $`0.040\pm 0.0005,`$ (63)
and for the ‘reduced’ ones,

$$\kappa ^{}=0.97\pm 0.13,\kappa ^{\text{FN}}=0.15\pm 0.01.$$
(64)
All these do not differ significantly from the values (5153) and (56) obtained by using the unrestricted data set. In particular, the sign of $`\lambda ^{}`$ remains positive. We therefore conclude that, also at large distances, the minimally modified FN action (41) fails to describe the Yang–Mills ensemble of $`𝒏`$–fields.

## 6 Summary and discussion

We have performed a lattice test of the FN conjecture stating that low–energy $`SU(2)`$ Yang–Mills theory is equivalent to a Skyrme–type sigma model. More specifically, FN suggest that the knot solitons of their model might be related to the Yang–Mills glueball spectrum.

Using standard Monte Carlo techniques, we have generated an ensemble of $`SU(2)`$ link fields from the Wilson action. This ensemble was then used to extract an associated ensemble of color vectors $`𝒏`$, with $`𝒏`$ parametrizing the gauge invariant distance between the maximally–Abelian and lattice Landau gauge slices. As these gauges are close to each other, there is a preferred direction for the $`𝒏`$–field which corresponds to explicit symmetry breaking. In this way we avoid the appearance of massless Goldstone bosons and thus generate a nonvanishing mass gap. A study of the exponential decay of correlators yields a mass gap close to 1 GeV. To find the effective action describing the Yang–Mills ensemble of $`𝒏`$–fields we have employed inverse Monte Carlo techniques. These are based on Schwinger–Dyson and Ward identities which we have derived analytically on the lattice. The identities have been evaluated numerically for the Yang–Mills ensemble on the one hand, and for ensembles stemming from LO and NLO effective actions on the other hand. As a result, we have found strong evidence that the ensemble generated from Yang–Mills theory cannot be described by the FN action plus a minimal symmetry–breaking term to allow for a mass gap. This follows from a number of discrepancies between the two ensembles. First, and most prominent, the sign of $`\lambda ^{}`$ is positive, implying negative magnetization $`𝔐`$, at variance with the value from the Yang–Mills ensemble. Second, the reduced two–point function $`h`$ ($`h^{\text{FN}}`$) from the NLO ensemble increases (decreases) with lattice disctance $`x`$, while for Yang–Mills the behavior is just the opposite. Third, the size of the mass gap is larger than for the Yang–Mills ensemble of $`𝒏`$–fields.

It is quite conceivable that magnetization (and susceptibility) can be recovered correctly by adding more (symmetry–breaking) terms to the NLO action (work in this direction is under way). The same remark applies to the mass gap. Note, however, that one cannot naturally expect the Yang–Mills and $`\sigma `$–model mass gaps to coincide due to the nonlocal relation between $`𝒏`$ and the link variables $`U`$333 As there is no unique or natural definition for $`𝒏`$, one may try alternative prescriptions for $`𝒏=𝒏[U]`$. A fairly local one is the following. Write the (gauge fixed) links as $`U_{x,\mu }=u_{x,\mu }^0+iu_{x,\mu }`$. Then define $`\overline{n}_x\overline{u}_x/|\overline{u}_x|`$ with the link average $`\overline{u}_x_\mu u_{x,\mu }`$. Under global gauge transformations this transforms properly such that $`\overline{𝒏}`$ is another color unit vector.. Whether this represents a problem is a question of scales. If the effective $`\sigma `$–model were valid only for distances of, say, $`R`$ $`\stackrel{>}{}`$ 0.8 fm corresponding to energies $`E`$ $`\stackrel{<}{}`$ 0.25 GeV, as suggested by the discussion of Section 5, then it would make no sense to address questions like the glueball spectrum. An analogous situation holds for the Fermi theory of weak interactions which also is only effective much below the $`W`$ and $`Z`$ scales.

Finally, one should mention that there is still another fundamental problem associated with describing Yang–Mills theory in terms of effective $`\sigma `$–models. Allowing for finite temperature, the latter are in the universality class of the 4$`d`$ Heisenberg model, while $`SU(2)`$ Yang–Mills theory is in the 3$`d`$ Ising class . This issue has been discussed recently in the context of constructing effective actions via Abelian projections . Again, if the $`\sigma `$–model scale were below the critical temperature, the effective theory would only be valid in the confined phase and would have nothing to say about the behavior of Yang–Mills theory close to the phase transition. Otherwise, an infinite number of operators would be required which, of course, is anything else but an ‘effective’ description. Summarizing, we conclude that, while a reasonable effective model generalizing the FN action may exist in principle, it will be of little practical use.

###### Acknowledgments.

The authors are indebted to S. Shabanov for suggesting this investigation and to P. van Baal for raising the issue of Goldstone bosons. Discussions with F. Bruckmann, P. de Forcrand, A. di Giacomo, A. González-Arroyo, D. Hansen, M. Hasenbusch, K. Langfeld, E. Seiler, M. Teper, and V. Zakharov are gratefully acknowledged. The work of T.H. was supported by DFG under contract Wi 777/5-1. He thanks the theory department of LMU Munich, in particular V. Mukhanov, for the hospitality extended to him. L.D. thanks G. Bali for providing the code for algorithm AI as well as helpful advice.

## Appendix A Conventions

Left and right lattice derivatives are defined as

$`_\mu f_x`$ $``$ $`f_{x+\mu }f_x,`$ (65)
$`_\mu ^{}f_x`$ $``$ $`f_{x\mu }f_x.`$ (66)
The ordinary lattice Laplacian $`\mathrm{}_\mu ^{}_\mu `$ is a negative semi–definite operator. Its action on lattice functions $`f`$ is given by

$$\mathrm{}f_x\underset{\mu }{}(2f_xf_{x+\mu }f_{x\mu }).$$
(67)
The covariant Laplacian $`\mathrm{}[U]`$ in the adjoint representation acts as

$$\mathrm{}^{ab}[U]f_x^b\underset{\mu }{}\left(2f_x^aR_{x,\mu }^{ab}f_{x+\mu }^bR_{x\mu ,\mu }^{ba}f_{x\mu }^b\right),$$
(68)
where we have defined the adjoint link

$$R_{x,\mu }^{ab}\frac{1}{2}\text{tr}(\tau ^aU_{x,\mu }\tau ^bU_{x,\mu }^{}).$$
(69)
## Appendix B Relating LLG and MAG

From (9) it follows immediately that the LLG minimizes the functional

$$\overline{F}_{\mathrm{LLG}}\underset{l}{}\text{tr}(𝟙{}_{}{}^{\mathrm{\Omega }}U_{l}^{}),l(x,\mu ),$$
(70)
and thus tends to bring the links $`U_l`$ close to $`𝟙`$. The MAG, on the other hand, minimizes

$$\overline{F}_{\mathrm{MAG}}\underset{l}{}(1{}_{}{}^{g}R_{l}^{33}),$$
(71)
and thus wants to bring the 33–entry of the adjoint link $`R_l^{ab}`$ close to 1. From (69) it is obvious that, if $`U_l`$ equals unity, the same will be true for $`R_l`$. This can be made more precise: if $`U_l𝟙+iaA_l`$, $`A_l`$ hermitean, then it is an easy exercise to show that $`R_l^{33}=1+O(a^2)`$. In this sense, the LLG is close to the MAG.

## Appendix C Schwinger–Dyson equations and Ward identities

We begin with computing the infinitesimal rotations of the various contributions in (15) and (16) to the effective action. It turns out that, for all $`S_j`$, the action of the angular momentum can be written as

$$i𝑳_xS_j=𝒏_x\times 𝑲_{jx},$$
(72)
(and analogous for the $`S_k^{}`$) with the vectors $`𝑲_{jx}`$ and $`𝑲_{kx}^{}`$ given by

$`𝑲_{1x}`$ $`=`$ $`2\mathrm{}𝒏_x`$ (73)
$`𝑲_{2x}`$ $`=`$ $`2\mathrm{}^2𝒏_x`$ (74)
$`𝑲_{3x}`$ $`=`$ $`2\left[\mathrm{}𝒏_x(𝒏_x\mathrm{}𝒏_x)+\mathrm{}\left(𝒏_x(𝒏_x\mathrm{}𝒏_x)\right)\right]`$ (75)
$`𝑲_{4x}`$ $`=`$ $`2\left[_\mu ^{}_\nu 𝒏_x(𝒏_x_\mu ^{}_\nu 𝒏_x)+_\mu ^{}_\nu \left(𝒏_x(𝒏_x_\mu ^{}_\nu 𝒏_x)\right)\right]`$ (76)
$`𝑲_{1x}^{}`$ $`=`$ $`𝒉`$ (77)
$`𝑲_{2x}^{}`$ $`=`$ $`2𝒉(𝒏_x𝒉)`$ (78)
$`𝑲_{3x}^{}`$ $`=`$ $`𝒉(𝒏_x\mathrm{}𝒏_x)+\mathrm{}𝒏_x(𝒏_x𝒉)+\mathrm{}(𝒏_x(𝒏_x𝒉)).`$ (79)
Choosing the $`F`$’s in (26) as $`n_x^a`$, $`n_x^an_y^b`$ and $`n_x^an_y^bn_z^c`$, respectively, results in the Ward identities

$`G_{\underset{¯}{x}y}^{ii}\lambda _1^{}+2G_{\underset{¯}{x}\underset{¯}{x}y}^{3ii}\lambda _2^{}+G_{\underset{¯}{x}\underset{¯}{x}y}^{\mathrm{\Delta }ii}\lambda _3^{}`$ $`=`$ $`2G_y^3,`$ (80)
$`G_{\underset{¯}{x}yz}^{ii3}\lambda _1^{}+2G_{\underset{¯}{x}\underset{¯}{x}yz}^{3ii3}\lambda _2^{}+G_{\underset{¯}{x}\underset{¯}{x}yz}^{\mathrm{\Delta }ii3}\lambda _3^{}`$ $`=`$ $`2G_{yz}^{33}+G_{yz}^{ii},`$ (81)
$`G_{\underset{¯}{x}yzz^{}}^{ii33}\lambda _1^{}+2G_{\underset{¯}{x}\underset{¯}{x}yzz^{}}^{3ii33}\lambda _2^{}+G_{\underset{¯}{x}\underset{¯}{x}yzz^{}}^{\mathrm{\Delta }ii33}\lambda _3^{}`$ $`=`$ $`2G_{yzz^{}}^{333}+2G_{yzz^{}}^{i(i3)}.`$ (82)
Here, the superscript $`(i3)`$ denotes symmetrization in $`i`$, 3, and we have introduced the shorthand notations

$`G_{xyz\mathrm{}}^{abc\mathrm{}}`$ $``$ $`n_x^an_y^bn_z^c\mathrm{},`$ (83)
$`G_{\underset{¯}{x}y\mathrm{}}^{ii\mathrm{}}`$ $``$ $`{\displaystyle \underset{x}{}}G_{xy\mathrm{}}^{ii\mathrm{}},`$ (84)
$`G_{\underset{¯}{x}\underset{¯}{x}y\mathrm{}}^{3ii\mathrm{}}`$ $``$ $`{\displaystyle \underset{x}{}}n_x^3n_x^in_y^i\mathrm{},`$ (85)
$`G_{\underset{¯}{x}\underset{¯}{x}y\mathrm{}}^{\mathrm{\Delta }ii\mathrm{}}`$ $``$ $`{\displaystyle \underset{x}{}}(n_x^a\mathrm{}n_x^a)n_x^in_y^i\mathrm{},`$ (86)
with $`\underset{¯}{x}`$ denoting summation over all lattice sites $`x`$. In particular, one has

$$G_y^3𝔐,G_{\underset{¯}{x}y}^{ii}2\chi ^{},G_{xy}^{ii}2G_{xy}^{},G_{xy}^{33}G_{xy}^{}.$$
(87)
It seems that a particular recurrence pattern arises in (8082) that could be used to derive a Ward identity for an insertion of an arbitrary number of $`𝒏`$’s. For the NLO derivative expansion, however, the three identities are sufficient to determine the symmetry–breaking couplings $`\lambda ^{}`$.

For the particular case that $`F`$ in (23) equals $`n`$, $`F_x^a=n_x^a`$, we can give the general Schwinger–Dyson equation in closed form,

$$\underset{j}{}\lambda _jH_{j,xy}+\underset{k}{}\lambda _k^{}H_{k,xy}^{}=𝔐\delta _{xy},$$
(88)
where we have introduced the two–point function

$$H_{j,xy}n_x^in_y^{[i}K_{jy}^{3]},$$
(89)
and, completely analogous, $`H_{k,xy}^{}`$. The FN terms are (75) minus (76), so that the relevant two–point function becomes

$$H_{xy}^{\text{FN}}H_{3,xy}H_{4,xy},$$
(90)
which has been used in (42).