text:
# Quantum Frenkel-Kontorova Model

## 1 Introduction

Many physical systems such as charge density waves, magnetic spirals, Josephson junctions and adsorbed monolayers exhibit a competition between two or more length scales. To understand such modulated structures, a simple model, known as the Frenkel-Kontorova (FK) model , was introduced more than half a century ago. The FK model describes a chain of atoms connected by springs in the presence of an external sinusoidal potential. It has two length scales: the natural length of the spring and the period of the external potential. Due to the competition between these two length scales, the FK model exhibits a wealth of interesting and complex phenomena. Extensive studies of the FK model have been made since its introduction. In most of the earlier studies, the FK model was treated in the continuum approximation. Although the continuum approximation, which reproduces the sine-Gordon equation and its soliton solutions, gives some very useful information, it misses many essential features. A major advance in the study of the FK model was made by Aubry . Aubry reverted to the discrete version of the FK model and made essential use of the seminal Kolmogorov-Arnol’d-Moser (KAM) theorem. In analogy to the breakup of a KAM torus, Aubry discovered the “transition by breaking of analyticity.” The Aubry transition is similar to a phase transition, thus all the apparatus invented in the study of second-order phase transitions can be borrowed to study this transition. A fairly complete picture of the FK model was thus obtained.

So far most of the studies made on the FK model have been confined to the classical model. Very little study has been done on the quantum FK model. Although the FK model looks simple, it’s nevertheless a non-integrable interacting many-body system. There aren’t that many methods at our disposal to study the quantum behavior of such a model. Moreover, besides theoretical interest, there are also practical considerations that demand a deeper understanding of the quantum FK model. For example, the study of friction, or tribology , is a very important subject both in terms of fundamental science and applications. And indeed, the FK model has been used to study tribology. However, as technology is moving into the nano-regime, nanotribology  is also becoming a very important subject of study. In the nano-regime, it’s more likely than not that quantum mechanics will play an important role. This makes an understanding of the quantum FK model not only desirable but very likely imperative. This is the motivation for our interest in the quantum FK model.

The first study of the quantum FK model was made by Borgonovi, Guarneri and Shepelyansky . They mainly used a numerical approach. Later, Berman, Bulgakov and Campbell  studied the same problem using a mean-field approach. In our study we will use a hybrid analytical and numerical approach. The analytical method is based on the squeezed state theory, and the numerical method is based on the quantum Monte Carlo (QMC) algorithm.

This paper gives a review of our recent work on the quantum FK model . It is organized as follows. In Sec. 2, we will define the quantum Frenkel-Kontorova model and use the Monte-Carlo method to study it. In Sec. 3 we will formulate the squeezed state approach. In Sec. 4 we will study the quantum effects on the Aubry transition. In Sec. 5 some concluding remarks will be given.

## 2 Quantum Frenkel-Kontorova model

The Hamiltonian of the quantum Frenkel-Kontorova model is described by:

$$\widehat{}=\underset{i}{}\left[\frac{\widehat{p}_i^2}{2m}+\frac{\gamma }{2}(\widehat{x}_{i+1}\widehat{x}_ia)^2V\mathrm{cos}(q_0\widehat{x}_i)\right].$$
(1)
Here $`m`$ is the mass of the particles, $`\gamma `$ the elastic constant of the spring, $`2\pi /q_0`$ the period of the external potential, $`V`$ the magnitude of the external potential, and $`a`$ the equilibrium distance between two nearest neighbor particles. For convenience, we shall rescale the variables so as to obtain a dimensionless Hamiltonian

$$\widehat{H}=\underset{i}{}\left[\frac{\widehat{P}_i^2}{2}+\frac{1}{2}(\widehat{X}_{i+1}\widehat{X}_i\mu )^2K\mathrm{cos}\widehat{X}_i\right],$$
(2)
where $`\widehat{X}_i=q_0\widehat{x}_i`$, $`\widehat{P}_i=(\frac{q_0}{\sqrt{m\gamma }})\widehat{p}_i`$, $`\mu =q_0a`$, and $`K=Vq_0^2/\gamma `$. We define an “effective Planck’s constant”:

$$\stackrel{~}{\mathrm{}}=\left(\frac{q_0^2}{\sqrt{m\gamma }}\right)\mathrm{}.$$
(3)
This effective Planck’s constant is the ratio of the natural quantum energy scale $`\mathrm{}\omega _0`$ to the natural classical energy scale $`\gamma /q_0^2`$, and $`\omega _0^2=\gamma /m`$. This quantity provides a measure of the importance of quantum effects. In particular, quantum effects become crucial at zero temperature when thermal fluctuations vanish.

In the classical FK model, the positions of the particles in the ground state, denoted by $`X_i^e`$, are determined by the following equation:

$$X_{i+1}^e2X_i^e+X_{i1}^e=K\mathrm{sin}X_i^e.$$
(4)
This equation can be cast into the form of the standard map by defining $`I_{i+1}=X_{i+1}^eX_i^e`$:

$$\begin{array}{c}I_{i+1}=I_i+K\mathrm{sin}X_i^e,\hfill \\ X_{i+1}^e=I_{i+1}+X_i^e.\hfill \end{array}$$
(5)
To describe the transition from the unpinned phase in the subcritical regime to the pinned phase in the supercritical regime, the hull function is introduced. It is defined by

$$X_i^e=f(X_i).$$
(6)
Here $`X_i`$ is the unperturbed ($`K=0`$) particle position. Another useful function is called the $`g`$-function, defined by

$$g=\frac{1}{K}(X_{i+1}^e2X_i^e+X_{i1}^e)$$
(7)
For $`K<K_c`$, the ground state of the classical FK model corresponds to an invariant curve in the standard map (5). The hull function is a monotonic analytical function \[Fig. 1(a)\], and the $`g`$-function is $`g(X)=\mathrm{sin}X`$ \[Fig. 1(b)\]. At $`K=K_c`$, the last invariant curve in the standard map breaks up and becomes a Cantorus. The hull functions becomes a monotonic function with a countable number of discontinuities \[Fig. 1(c)\], and the $`g`$-function becomes a Cantor subset of the function $`\mathrm{sin}X`$ \[Fig. 1(d)\].

By using the quantum Monte Carlo method, we can calculate the ground state averages of the particle positions $`\overline{X}_i`$ and plot the following three functions so as to compare with their classical counterparts: (1) The quantum hull function, i.e., $`\overline{X}_i`$ as a function of the unperturbed positions $`X_i`$; (2) the quantum $`g`$-function which is defined as

$$g_q=\frac{1}{K}(\overline{X}_{i+1}2\overline{X}_i+\overline{X}_{i1});$$
(8)
and (3) ($`\overline{X}_i,\overline{I}_i`$) in the phase space of the standard map.

The quantum effect does not destroy the functional dependence of $`g_q`$ on $`\overline{X}_i`$, instead it renormalizes it to a sawtooth map \[Fig. 2\]. This phenomenon was first observed by Borgonovi et al.  and was interpreted as a consequence of quantum tunneling between the edges of the gap in the hull function.

## 3 Squeezed state approach

In this section, we will study the effect of quantum fluctuations on the FK model by using the squeezed state approach .

The squeezed state was first used in quantum optics and has proven to be very useful in many other fields . The squeezed state approach starts from the time-dependent variational principle

$$\delta 𝑑t\mathrm{\Psi }(t)|i\mathrm{}\frac{}{t}\widehat{H}|\mathrm{\Psi }(t)=0.$$
(9)
Variation with respect to (w.r.t.) $`\mathrm{\Psi }(t)|`$ gives rise to the Schrödinger equation. The exact solution may be approximated by restricting the choice of states to a subspace of the full Hilbert space. In the squeezed state approach, the squeezed state is chosen as $`|\mathrm{\Psi }(t)`$.

The squeezed state $`|\mathrm{\Phi }(\alpha ,\beta )`$ is defined by

$$|\mathrm{\Phi }(\alpha ,\beta )=e^{\widehat{S}(\alpha )}e^{\widehat{T}(\beta )}|0,$$
(10)
where

$$\widehat{S}(\alpha )=\underset{i}{}(\alpha _i\widehat{a}_i^{}\alpha _i^{}\widehat{a}_i)$$
(11)
is the Weyl’s displacement operator, and

$$\widehat{T}(\beta )=\frac{1}{2}\underset{ij}{}(\widehat{a}_i^{}\beta _{ij}\widehat{a}_j^{}\widehat{a}_i\beta _{ij}^{}\widehat{a}_j).$$
(12)
is the squeezed state operator. $`\widehat{S}^{}(\alpha )=\widehat{S}(\alpha ),\widehat{T}^{}(\beta )=\widehat{T}(\beta )`$. For brevitiy, we will use the abbreviation $`|\mathrm{\Phi }|\mathrm{\Phi }(\alpha ,\beta )`$. If $`\beta =0`$, the squeezed state is reduced to the coherent state. However, as we shall see later, the coherent state is not adequate to study quantum fluctuations.

The position and momentum operators for the $`i`$th particle are written as

$$\begin{array}{c}\widehat{X}_i=\sqrt{\frac{\stackrel{~}{\mathrm{}}}{2}}(\widehat{a}_i^{}+\widehat{a}_i),\hfill \\ \widehat{P}_i=i\sqrt{\frac{\stackrel{~}{\mathrm{}}}{2}}(\widehat{a}_i^{}\widehat{a}_i).\hfill \end{array}$$
(13)
Here, $`\widehat{a}_i^{}`$ and $`\widehat{a}_i`$ are the creation and annihilation operators which satisfy the canonical commutation relations: $`[\widehat{a}_i,\widehat{a}_j^{}]=\delta _{ij}`$, $`[\widehat{a}_i,\widehat{a}_j]=0`$ and $`[\widehat{a}_i^{},\widehat{a}_j^{}]=0`$.

Using $`|\mathrm{\Phi }`$ as a trial wave function, we can find the expectation values of the position and momentum operators of the $`i`$th particle :

$$\begin{array}{c}\overline{X}_i\mathrm{\Phi }|\widehat{X}_i|\mathrm{\Phi }=\sqrt{\frac{\stackrel{~}{\mathrm{}}}{2}}(\alpha _i^{}+\alpha _i),\hfill \\ \overline{P_i}\mathrm{\Phi }|\widehat{P}_i|\mathrm{\Phi }=i\sqrt{\frac{\stackrel{~}{\mathrm{}}}{2}}(\alpha _i^{}\alpha _i).\hfill \end{array}$$
(14)
Fluctuations in the position and the momentum are given by

$`\mathrm{\Delta }X_i^2`$ $``$ $`\mathrm{\Phi }|(\widehat{X}_i\overline{X}_i)^2|\mathrm{\Phi }=\stackrel{~}{\mathrm{}}G_{ii},`$
$`\mathrm{\Delta }P_i^2`$ $``$ $`\mathrm{\Phi }|(\widehat{P}_i\overline{P}_i)^2|\mathrm{\Phi }=\stackrel{~}{\mathrm{}}\left({\displaystyle \frac{G_{ii}^1}{4}}+4{\displaystyle \underset{l,k}{}}\mathrm{\Pi }_{il}G_{lk}\mathrm{\Pi }_{ki}\right).`$ (15)
The fluctuation covariance between the $`i`$th particle and the $`j`$th particle is

$$\mathrm{\Delta }X_i\mathrm{\Delta }X_j\mathrm{\Phi }|(\widehat{X}_i\overline{X}_i)(\widehat{X}_j\overline{X}_j)|\mathrm{\Phi }=\stackrel{~}{\mathrm{}}G_{ij},$$
(16)
where

$$\begin{array}{c}G_{ij}=\frac{1}{2}(\mathrm{cosh}^2\sqrt{\beta \beta ^{}}+\mathrm{sinh}^2\sqrt{\beta \beta ^{}})_{ij}+\frac{1}{2}(M\beta +\beta ^{}M)_{ij},\hfill \\ \mathrm{\Pi }_{ij}=\frac{i}{4}G_{ij}^1(M\beta \beta ^{}M)_{ij},\hfill \end{array}$$
(17)
and,

$$M=\frac{\mathrm{sinh}\sqrt{\beta \beta ^{}}\mathrm{cosh}\sqrt{\beta \beta ^{}}}{\sqrt{\beta \beta ^{}}}.$$
(18)
Since $`\beta `$ is a symmetric matrix, $`G_{ij}=G_{ji}`$ and $`\mathrm{\Pi }_{ij}=\mathrm{\Pi }_{ji}`$. Using the following relation

$$\mathrm{\Phi }|\mathrm{cos}\widehat{X}_i|\mathrm{\Phi }=\mathrm{exp}\left(\frac{\stackrel{~}{\mathrm{}}}{2}G_{ii}\right)\mathrm{cos}\overline{X}_i,$$
(19)
we obtain the expectation value of the Hamiltonian, $`\overline{H}\mathrm{\Phi }|\widehat{H}|\mathrm{\Phi }`$,

$`\overline{H}`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{1}{2}}\left(\overline{P}_i^2+\stackrel{~}{\mathrm{}}({\displaystyle \frac{G_{ii}^1}{4}}+4{\displaystyle \underset{l,k}{}}\mathrm{\Pi }_{il}G_{lk}\mathrm{\Pi }_{ki})\right)+{\displaystyle \underset{i}{}}{\displaystyle \frac{1}{2}}(\overline{X}_{i+1}\overline{X}_i\mu )^2`$ (20)
$`+`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{\stackrel{~}{\mathrm{}}}{2}}\left(G_{ii}+G_{i+1i+1}2G_{i+1i}\right){\displaystyle \underset{i}{}}K\mathrm{exp}\left({\displaystyle \frac{\stackrel{~}{\mathrm{}}}{2}}G_{ii}\right)\mathrm{cos}\overline{X}_i.`$
The variables $`\overline{X}_i`$ and $`\overline{P}_i`$, and $`G_{ij}`$ and $`\mathrm{\Pi }_{ij}`$ form explicitly canonical conjugates . To find the ground state of the quantum FK model, we will use a variational approach in which these four variables are treated as independent variables. Variation with respect to $`\overline{P}_i`$ and $`\overline{X}_i`$ yields, respectively,

$$\begin{array}{c}\overline{P}_i=0,\hfill \\ \overline{X}_{i+1}2\overline{X}_i+\overline{X}_{i1}=K_i\mathrm{sin}\overline{X}_i,\hfill \end{array}$$
(21)
where

$$K_i=K\mathrm{exp}\left(\frac{\stackrel{~}{\mathrm{}}}{2}G_{ii}\right).$$
(22)
This equation determines the expectation value of the particle’s position. Unlike its classical counterpart ($`\stackrel{~}{\mathrm{}}=0`$, $`K_i=K`$), this equation is coupled to the quantum fluctuation by $`\stackrel{~}{\mathrm{}}G_{ii}`$. Variation w.r.t. $`\mathrm{\Pi }_{ij}`$ leads to

$$\underset{k}{}G_{ik}\mathrm{\Pi }_{kj}=0.$$
(23)
To obtain the equation for $`G_{ij}`$, we first take variation w.r.t. $`G_{ik}`$ and note that

$$\frac{\delta G_{ij}}{\delta G_{kl}}=\delta _{ik}\delta _{jl}.$$
(24)
We then multiply both sides of the equation by $`G_{kj}`$ and sum over $`k`$. We then get the closed set of equations for the covariance $`G_{ij}`$:

$$(GF)_{ij}=G_{i1j}+G_{i+1j},$$
(25)
where

$$F_{ij}=\left(1+\frac{K_i}{2}\mathrm{cos}\overline{X}_i\right)\delta _{ij}\frac{(G^2)_{ij}}{8}.$$
(26)
This is a set of equations determining the quantum fluctuations of the particles. $`G`$ is an $`N\times N`$ symmetric matrix which provides information about fluctuations. Its diagonal elements give the variance of each particle while its off-diagonal elements give the covariance between particles. From this information we can calculate the correlation function of the quantum fluctuations. These equations are coupled to the expectation values $`\overline{X}_i`$.

Up to this point, we have obtained $`\frac{1}{2}N\times (N+1)+N`$ equations for all the variables. These equations provide some qualitative information about the system. If we define $`\overline{I}_{i+1}=\overline{X}_{i+1}\overline{X}_i`$, Eq. (21) can be cast into the form of a standard map,

$$\begin{array}{c}\overline{I}_{i+1}=\overline{I}_i+K_i\mathrm{sin}\overline{X}_i,\hfill \\ \overline{X}_{i+1}=\overline{I}_{i+1}+\overline{X}_i.\hfill \end{array}$$
(27)
Similarly, if we define $`Q_{i+1j}=G_{i+1j}G_{ij}`$, Eq.(25) can also be cast into the form of a map:

$$\begin{array}{c}Q_{i+1j}=Q_{ij}+[G(F2)]_{ij},\hfill \\ G_{i+1j}=G_{ij}+Q_{i+1j}.\hfill \end{array}$$
(28)
The difference between the classical case ($`\stackrel{~}{\mathrm{}}=0`$) and the quantum case is readily seen in Eq.(27) . In the classical case the coupling constant does not change from potential to potential; however, in the quantum case, the effective coupling constant changes from particle to particle. Because $`G_{ii}>0`$ for any $`\stackrel{~}{\mathrm{}}0`$, $`K_i<K`$, which implies that the quantum fluctuation will reduce the coupling strength. Another important difference is that in the classical case the positions of the particles in the ground state are determined by the standard map whereas in the quantum case they are determined by $`(N+1)`$-coupled two-dimensional maps. This makes the quantum FK model extremely difficult to study analytically.

One notices that when $`\beta =0`$, $`G_{ij}=1/2`$, which is the result of the coherent state theory. This is not the case in the quantum FK model. In Fig. 4 we plot $`\stackrel{~}{\mathrm{}}G_{ij}`$ as a function of $`i`$ and $`j`$. The data are obtained from QMC. It is obvious that $`G_{ij}1/2`$. Thus it seems that the coherent state is not adequate for the study of the quantum FK model.

We now make some comparisons with the QMC results. As mentioned before, to find the solution of these two sets of equations \[Eqs.(21) and (25)\] is equivalent to finding a periodic orbit in the $`2(N+1)`$-dimensional map \[Eqs.(27) and (28)\]. This is a very difficult problem to solve even numerically. Nevertheless, we will try to solve Eq.(21) self-consistently to see whether this equation can indeed give rise to the sawtooth map .

Using QMC, we obtained the expectation value of the particle’s position $`\overline{X}_i`$, from which we can construct the quantum hull function $`f_q`$, as shown in Fig. 2(a) for $`K=0.5`$ and (c) for $`K=5`$, respectively. The quantum $`g_q`$-function is shown in Fig. 2(b) and (d) for $`K=0.5`$ and $`K=5`$, respectively. The quantum fluctuation $`\stackrel{~}{\mathrm{}}G_{ii}`$ calculated from QMC is shown in Fig. 3 for the case $`K=5`$.

As usual, we use $`Q`$ particles and $`P`$ external potentials with period $`2\pi `$, and the periodic boundary condition is imposed : $`\overline{X}_{Q+i}=\overline{X}_i+2\pi P`$. The results shown in the figure are for $`P/Q=34/55`$.

To compare the squeezed state results with those from QMC, we substitute $`G_{ii}`$ calculated from QMC into Eq.(21). We then compute the expectation value of the particles’ positions by using Aubry’s gradient method. We then construct the quantum hull function $`f_q`$ and the quantum $`g_q`$-function shown in Fig. 5. The results agree surprisingly well with those from QMC. The most striking feature is the sawtooth shape of the $`g_q`$-function. This phenomenon was first observed by Borgonovi et al.  in their QMC computation and has been explained as a tunneling effect. Later, Berman et al.  obtained this result by using a mean field theory. In the framework of the squeezed state theory, the renormalization of the standard map to a sawtooth map is a consequence of quantum fluctuations. Our results show that the squeezed state approach indeed captures the effects of quantum fluctuations.

## 4 Quantum effects on the Aubry transition

In this section we will study the effects of quantum fluctuations on the Aubry transition. We will first study the behavior of the wave function of an incommensurate ground state. We will then introduce various indicators that might be used to characterize the quantum Aubry transition.

### 4.1 Ground state wave function of the quantum Frenkel-Kontorova model

We will use the QMC  to study the ground state of Hamiltonian (2). As in the classical case, we will concentrate our attention on the incommensurate state characterized by the golden mean winding number $`\omega _G=(\sqrt{5}1)/2`$. In the classical case the winding number is defined as the average separation of atoms per period, i.e., $`\omega =lim_N\mathrm{}\frac{X_NX_0}{2N\pi }`$. If $`\omega `$ is a rational number, the ground state is commensurate; and if it is an irrational number, the ground state is incommensurate. In the quantum case, $`\overline{X}_i`$ is defined as the expectation value of the position of the $`i`$th particle. As usual, we use the method of continued fraction expansion and approximate $`\omega _G`$ by its rational convergents $`\omega _n=F_n/F_{n+1}`$, where $`F_n`$ is the $`n`$th Fibonacci number defined by the recursion relation $`F_{n+1}=F_n+F_{n1}`$ with $`F_0=0,F_1=1`$. In our computation, we choose $`F_{n+1}`$ particles embedded in $`F_n`$ external potentials.

Since the external potential is periodic with period $`2\pi `$, we can fold the wave function to this period and then take the average over all the particles in the interval $`[0,2\pi ]`$. This quantity gives the probability of finding the particle at a given position $`X`$. We plot the averaged probability density $`|\mathrm{\Psi }|^2`$ in Fig. 6 with 144 particles embedded in 89 potentials for a fixed $`\stackrel{~}{\mathrm{}}`$ ($`=0.2`$) but different values of $`K`$. Here $`|\mathrm{\Psi }|^2`$ is normalized, i.e., $`_0^{2\pi }𝑑X|\mathrm{\Psi }(X)|^2=1`$. We observe that, in the small $`K`$ regime, the probability of finding the particles at any point of the potential is almost the same (see, for example, the curve for $`K=0.1`$). This is quite similar to the unpinned phase in the classical case. However, as the coupling constant increases, the probability of finding the particle in the upper part of the potential decreases whereas the probability of finding the particle in the lower part of the potential increases, as is noticeable from the appearance of peaks in the curves. As the coupling constant increases, the probability of finding the particles at the top is almost zero (see, for example, the curves for $`K=2`$ and $`5`$).

### 4.2 Indicators of the Aubry transition

#### 4.2.1 Disorder Parameter

In the classical FK model, Coppersmith and Fisher  have proposed a “disorder parameter” to describe the transition from the pinned phase to the unpinned phase. This parameter is defined as the minimum distance of a particle from the top of the potential, $`D=\mathrm{min}_{i,n}|X_i^e2\pi (n+\frac{1}{2})|`$. If $`D0`$, the particles are pinned; when $`D=0`$, they are unpinned. This “disorder parameter” also measures the discontinuity (or width of the biggest gap) of the hull function. One might want to use the same function to describe the quantum crossover. For instance, one may define a very similar quantum disorder parameter $`D_q=\mathrm{min}_{i,n}|\overline{X}_i2\pi (n+\frac{1}{2})|`$, where $`\overline{X}_i`$ is the expectation value of the $`i`$th particle’s position. However, this parameter $`D_q`$ could not capture the crossover of the ground state wave function. The reason is that, in the quantum case, the particle can tunnel from one side of the potential to the other. Thus the gap in the classical hull function does not survive quantum fluctuations (see Refs.  and ). It turns out that for large $`K`$, $`D_q`$ might still be close to (or equal to) zero. Therefore, a new parameter is needed to describe the wave function crossover. To this end, we define the probability of finding the particles at the potential top as a quantum “disorder parameter”,

$$P_t=\frac{1}{N}\underset{n=0}{\overset{N1}{}}|\mathrm{\Psi }(X)|^2\delta \left[X2\pi (n+\frac{1}{2})\right]𝑑X.$$
(29)
In Figs. 7(a)-(b), we plot $`P_t`$ as a function of $`K`$ for different temperatures $`T_e`$ ($`T_e=\stackrel{~}{\mathrm{}}/\tau `$, where $`\tau `$ is the length of the imaginary time axis) and different system sizes. In Fig. 7(a), we fix $`\stackrel{~}{\mathrm{}}(=0.2)`$, $`\omega =13/21`$ and set the temperature $`T_e=0.2/30,0.2/120,`$, $`0.2/480`$. The convergence is quite fast as is seen from the figure. The sharp decrease of $`P_t`$ is very evident for $`1<K<2`$. In this small $`K`$ regime, $`P_t`$ changes significantly: it drops almost two orders of magnitude. This drastic change can also be seen from Fig. 7(b). There the temperature is fixed at $`T_e=0.2/120`$, which can be regarded as effectively zero temperature. The system size is changed from 21 particles to 144 particles.

#### 4.2.2 Maximum Variance

Another parameter that can be used to depict this crossover is the maximal fluctuation of the particles. In the quantum case, we have observed that the particle situated at a position nearest to the top of the potential always has a maximal fluctuation since it has a higher classical potential energy. This observation has been verified numerically by QMC and theoretically by the squeezed state theory . Thus we can use this maximal fluctuation as another measure of quantum “disorder”. To make a qualitative comparison with the classical disorder parameter $`D`$, we define

$$\mathrm{\Delta }=\underset{i}{\mathrm{max}}\left[\sqrt{(X_i\overline{X}_i)^2}\right].$$
(30)
Here the average $`\mathrm{}`$ is taken over all the paths (approximately 4,000) used in QMC. $`\mathrm{\Delta }`$ as a function of $`K`$ is plotted in Fig. 8. The computations given in this figure have been carried out with $`\omega _n=F_n/F_{n+1}=34/55`$.

It is noticeable that the transition of the wave function is characterized by the different $`K`$-dependent behavior of $`\mathrm{\Delta }`$. In the small $`K`$ regime, $`\mathrm{\Delta }`$ is a constant independent of $`K`$ but changes with $`\stackrel{~}{\mathrm{}}`$. In the large $`K`$ regime, $`\mathrm{\Delta }`$ increases with $`K`$ but does not change with $`\stackrel{~}{\mathrm{}}`$. Furthermore, $`\mathrm{\Delta }`$ is approximately equal to the classical disorder parameter $`D`$. For comparison, we have included $`D`$ in the inset of Figure 8. In the classical case, the transition of course occurs at $`K=K_c`$.

These results show that the width and the shape of the probability density do not change with $`K`$ in the small $`K`$ regime. The width only depends on the strength of the quantum fluctuation $`\stackrel{~}{\mathrm{}}`$. However, this picture changes significantly for large $`K`$. In this regime, the profile of the probability density spreads out, and the width of the probability density is quite insensitive to quantum fluctuations and depends only on the coupling constant. It should however be stressed that the analogy between $`D`$ and $`\mathrm{\Delta }`$ cannot be carried too far since $`\mathrm{\Delta }`$ is not exactly the quantum correspondence of $`D`$. This is why even if we let $`\stackrel{~}{\mathrm{}}`$ go to zero, $`\mathrm{\Delta }`$ would not be zero in the small $`K`$ regime.

#### 4.2.3 Correlation function

As another indication of the quantum signature of the Aubry transition, we study the correlation function between the particles’ fluctuations. The cross-correlation function $`C_{ij}`$ is defined by

$$C_{ij}=\frac{l_{ij}}{\sqrt{l_{ii}l_{jj}}}.$$
(31)
It describes the fluctuation correlation between the particles $`i`$ and $`j`$. $`l_{ij}=(X_i\overline{X}_i)(X_j\overline{X}_j)`$ is the covariance of the fluctuations between the $`i`$th particle and the $`j`$th particle. Here $`\mathrm{}`$ denotes the ensemble average over space and ‘time’ in QMC . Since $`l_{ij}`$ is symmetric, $`C_{ij}=C_{ji}`$. In our numerical calculations, we take 4000 paths and each path is divided into many steps in time, ranging from 1500 to 3000 depending on $`\stackrel{~}{\mathrm{}}`$. The convergence is guaranteed by increasing the number of steps and the number of paths until the same results are reached.

Figure 9 shows the contour plots of $`C_{ij}`$ for $`K=0.5`$ and $`K=5`$. The plots are drawn in 12 equal steps from $`0.15`$ to $`1`$. A negative correlation function implies that the particles’ quantum fluctuations are out of phase. The sharp change of $`C_{ij}`$ from $`K=0.5`$ to $`K=5`$ is very readily noted. In the large $`K`$ regime, the correlation function matrix $`C_{ij}`$ is almost diagonal with a very small band width whereas in the small $`K`$ regime the number of nonzero off-diagonal elements increases and the band width is much larger. Thus, in the small $`K`$ regime, we have a much longer range interaction than in the large $`K`$ regime. This can also be seen from the averaged correlation function $`C(|ij|)=C_{ij}`$, where $`\mathrm{}`$ denotes the ensemble average over particles. We have investigated this function for a wide range of $`K`$. We tried to use different forms to best-fit this function, but found that it is very hard to see whether the decay is exponential, power-law or logarithmic. It seems that it obeys different laws in different ranges. Nevertheless, the correlation function in the large $`K`$ regime decays faster than that in the small $`K`$ regime, as illustrated in Fig. 10 for $`K=0.5`$ and $`K=5`$, respectively.

On the other hand, due to discreteness, it is also very difficult to define a parameter such as the width of the correlation function. As an alternative, we chose $`C(1)`$, which describes the correlation between the nearest-neighbor particles. In Fig. 11, we plot $`C^1(1)`$ against $`K`$ on a logarithmic scale. The signature of the transition is very similar to Fig. 8. In the small $`K`$ regime, $`C(1)`$ is a constant (about 0.5); and in the large $`K`$ regime, it decreases with $`K`$. Since the transition is largely smeared by quantum fluctuations, one can’t define a critical $`K_c`$ as in the classical case.

We know that if the correlation function decays exponentially in the short range, then $`\mathrm{ln}C^1(1)=l^1`$, which is just the Lyapunov exponent $`\sigma `$ in the classical case , $`l`$ being the correlation length. Thus, the parameter $`C(1)`$ captures the short-range correlation. Fig. 11 mimics the classical plot of the correlation length proposed by Aubry . In the classical case, for $`K<K_c`$, $`l=\mathrm{}`$ and $`\sigma =0`$, whereas for $`K>K_c`$, $`l`$ is finite and $`\sigma >0`$. However, in the quantum FK model, we do not have infinitely long interaction in the small $`K`$ regime. Instead, we have a finite interaction range, but it is much longer than that in the large $`K`$ regime.

## 5 Concluding remarks

The Frenkel-Kontorova model is a simple model. However, despite its deceptively simple appearance, it exhibits very complex behavior and contains very rich physics. It has found applications in a wide range of fields, from nonlinear dynamics to condensed matter physics to tribology. However, most of the studies have so far been confined to the classical Frenkel-Kontorova model. Very little work has been done on the quantum Frenkel-Kontorova model. We think this is indeed quite regrettable. The quantum Frenkel-Kontorova model is not only interesting from a theoretical point of view but is also important from a practical point of view. The work reported in this review represents only a first step in this direction. We however hope that it will serve as an impetus to further studies. It is our conviction that a thorough understanding of the quantum Frenkel-Kontorova model will lead to important breakthroughs both in theory and applications.

We would like to thank S. Aubry, F. Borgonovi, D. K. Campbell, H. Chen, R. B. Griffiths, C.-L. Ho, H.-Q. Lin, D. Shepelyansky, L.-H. Tang, and W.-M. Zhang for many helpful discussions. This work was supported in part by grants from the Hong Kong Research Grants Council (RGC) and the Hong Kong Baptist University Faculty Research Grant (FRG).


format: text

text:
# The evolution and persistence of dumbbells

## 1 Introduction

The formation of topological defects in any phase transition depends on the existence of a non-trivial, low-order homotopy group of the broken-symmetry vacuum manifold . Even in the absence of this, dynamically stable non-topological defects are sometimes still possible, but are usually assumed to be too weakly stable to lead to a lasting network. A particularly interesting case is non-topological strings where at least one example is known (semilocal strings ) in which a space-spanning network of strings can develop from the growth and joining of short string segments formed at the phase transition . Semilocal strings are a special case of electroweak strings .

The stability of electroweak strings has been analyzed in some detail in the Glashow-Salam-Weinberg (GSW) model. In the absence of fermions the theory is determined, up to scalings, by only two parameters: the weak mixing angle $`\theta _\mathrm{W}`$, and $`\beta `$ which is the square of the ratio of the Higgs and $`Z`$-boson masses (the observed value for the actual electroweak model is $`\mathrm{sin}^2\theta _\mathrm{W}0.23`$, and the precise value of $`\beta `$ has not yet been determined but is almost certainly greater than one).

In the case of interest here, the topology of the vacuum manifold is the three-sphere $`S^3`$, which does not support persistent topological defects in 3+1 dimensional spacetime, but there are grounds for expecting non-topological defects to form. In the limit $`\theta _\mathrm{W}=\pi /2`$ we recover the so-called semilocal model, whose $`S^3`$ vacuum manifold is known to support stable non-topological strings from both analytical  and numerical  work. In particular, a straight infinite semilocal string can be seen as a Nielsen-Olesen (NO) $`U(1)`$ vortex , embedded in a higher group $`SU(2)_{\mathrm{global}}\times U(1)_{\mathrm{local}}`$. Note, however, that in the semilocal string case, we can define a quantity (the winding number) which, although not a topological invariant in the usual sense, is topologically conserved.

The embedding of NO vortices into the full electroweak symmetry $`SU(2)_L\times U(1)_\mathrm{Y}`$ leads to $`Z`$-strings and $`W`$-strings (see ref. ). Although $`W`$-strings are expected to be unstable, analysis of infinite axially-symmetric $`Z`$-strings has shown the existence of a parameter regime where they are perturbatively stable . But these $`Z`$-strings are genuinely non-topological, there is no quantity which is topologically conserved.

In this paper we will consider the evolution and persistence of these genuine non-topological defects in a generalized GSW model, spanning all values of $`\beta `$ and $`\theta _\mathrm{W}`$, using numerical methods. Isolated, infinite, axially-symmetric $`Z`$-string configurations will not be formed in a realistic system, but configurations consisting of monopole-antimonopole pairs joined by $`Z`$-string segments, named *dumbbells* by Nambu , are certainly possible. Isolated dumbbells are expected to collapse under the string tension, at least in the absence of rotation, jittering or magnetic fields . Of particular interest is the question of whether there is a region of parameter space in which the density of dumbbells is sufficiently high so that they are able to generate a persistent network of strings by building up longer segments from the merging of shorter ones due to interaction of neighbouring monopoles, as has been found in the semilocal case . This region of parameter space will lie far from the measured values for $`\theta _\mathrm{W}`$ and $`\beta `$, but the mere existence of this region is striking. It shows that in models close to real physical ones, completely non-topological defect networks can form and persist, so in extensions of the GSW model (to higher symmetry groups or extra fields), or models with a background magnetic field , or even models in which topological (semilocal) defects exist in some limit, non-topological defects cannot be ruled out immediately.

Moreover, the study of the dynamics of dumbbells can be of interest in computing other early universe features, such as the primordial magnetic field helicity .

## 2 The model

The bosonic sector of the GSW electroweak model describes an $`SU(2)_L\times U(1)_\mathrm{Y}`$ invariant theory with a scalar field $`\mathrm{\Phi }`$ in the fundamental representation of $`SU(2)_L`$, with lagrangian

$$=\left|D_\mu \mathrm{\Phi }\right|^2\frac{1}{4}W_{\mu \nu }^aW^{a\mu \nu }\frac{1}{4}Y_{\mu \nu }Y^{\mu \nu }\lambda \left(\mathrm{\Phi }\mathrm{\Phi }^{}\frac{\eta ^2}{2}\right)^2.$$
(1)
The covariant derivative is given by

$$D_\mu _\mu \frac{ig_\mathrm{W}}{2}\tau ^aW_\mu ^a\frac{ig_\mathrm{Y}}{2}Y_\mu ,a=1,2,3,$$
(2)
where $`\mathrm{\Phi }`$ is a complex doublet, $`\tau ^a`$ are the Pauli matrices, $`W_\mu ^a`$ is a $`SU(2)`$ gauge field and $`Y_\mu `$ is a $`U(1)`$ gauge field. The field strengths associated with these gauge fields are

$`W_{\mu \nu }^a`$ $`=`$ $`_\mu W_\nu ^a_\nu W_\mu ^a+g_\mathrm{W}ϵ^{abc}W_\mu ^bW_\nu ^c;`$
$`Y_{\mu \nu }`$ $`=`$ $`_\mu Y_\nu _\nu Y_\mu ,`$ (3)
respectively, and there is no distinction between upper and lower group indices ($`ϵ^{123}=1`$).

When the scalar field acquires a non-zero vacuum expectation value the symmetry breaks from $`SU(2)_L\times U(1)_\mathrm{Y}`$ to $`U(1)_{\mathrm{e}.\mathrm{m}.}`$, leaving a massive scalar field ($`m_H=\sqrt{2\lambda }\eta `$), a massless neutral photon ($`A_\mu `$), a massive neutral $`Z`$-boson ($`Z_\mu `$, $`m_\mathrm{z}=g_\mathrm{z}\eta /2l_v^1`$, where $`g_\mathrm{z}=\sqrt{g_\mathrm{Y}^2+g_\mathrm{W}^2}`$), and two massive charged $`W`$-bosons ($`W_\mu ^\pm `$, $`m_W=g_\mathrm{W}\eta /2`$).

We make the following rescaling

$$\mathrm{\Phi }\frac{\eta }{\sqrt{2}}\mathrm{\Phi },x_\mu \frac{\sqrt{2}}{g_\mathrm{z}\eta }x_\mu ,g_\mathrm{Y}Y_\mu \frac{g_\mathrm{z}\eta }{\sqrt{2}}Y_\mu ,g_\mathrm{W}W_\mu ^a\frac{g_\mathrm{z}\eta }{\sqrt{2}}W_\mu ^a,$$
(4)
to choose $`l_v`$ as the unit of length, $`\eta `$ as the unit of energy, and the $`Z`$-charge of the scalar field ($`g_\mathrm{z}`$) as the unit of charge (up to numerical factors). This brings the classical field equations to the form

$`D^\mu D_\mu \mathrm{\Phi }+{\displaystyle \frac{2\lambda }{g_\mathrm{z}^2}}\left(\mathrm{\Phi }\mathrm{\Phi }^{}1\right)\mathrm{\Phi }`$ $`=`$ $`0;`$
$`_\nu W^{\mu \nu a}+ϵ^{abc}W_\nu ^bW^{\mu \nu c}`$ $`=`$ $`{\displaystyle \frac{i}{2}}\mathrm{cos}^2\theta _\mathrm{W}\left[\mathrm{\Phi }^{}\tau ^aD^\mu \mathrm{\Phi }\left(D^\mu \mathrm{\Phi }\right)^{}\tau ^a\mathrm{\Phi }\right];`$
$`_\nu Y^{\mu \nu }`$ $`=`$ $`{\displaystyle \frac{i}{2}}\mathrm{sin}^2\theta _\mathrm{W}\left[\mathrm{\Phi }^{}D^\mu \mathrm{\Phi }\left(D^\mu \mathrm{\Phi }\right)^{}\mathrm{\Phi }\right],`$ (5)
where the weak mixing angle is given by $`\mathrm{tan}\theta _\mathrm{W}g_\mathrm{Y}/g_\mathrm{W}`$, and now

$`W_{\mu \nu }^a`$ $``$ $`_\mu W_\nu ^a_\nu W_\mu ^a+ϵ^{abc}W_\mu ^bW_\nu ^c;`$
$`D_\mu `$ $``$ $`_\mu {\displaystyle \frac{i}{2}}\tau ^aW_\mu ^a{\displaystyle \frac{i}{2}}Y_\mu .`$ (6)
Typically the $`Z`$\- and $`A`$-fields are expressed in the unitary gauge $`\mathrm{\Phi }^T=(0,1)`$, but this choice is not well suited to work with defects. Instead, when there are points in space-time with $`|\mathrm{\Phi }|1`$, it is customary to use a more general definition of these fields which depends on the Higgs field configuration at each point , namely,

$`Z_\mu `$ $``$ $`\mathrm{cos}\theta _\mathrm{W}n^a(x)W_\mu ^a\mathrm{sin}\theta _\mathrm{W}Y_\mu ;`$
$`A_\mu `$ $``$ $`\mathrm{sin}\theta _\mathrm{W}n^a(x)W_\mu ^a+\mathrm{cos}\theta _\mathrm{W}Y_\mu ,`$ (7)
where

$$n^a(x)\frac{\mathrm{\Phi }^{}(x)\tau ^a\mathrm{\Phi }(x)}{\mathrm{\Phi }^{}(x)\mathrm{\Phi }(x)},$$
(8)
is a unit vector by virtue of the Fierz identity $`_a\left(\mathrm{\Phi }^{}\tau ^a\mathrm{\Phi }\right)^2=\left(\mathrm{\Phi }^{}\mathrm{\Phi }\right)^2`$. There are also several possible definitions for the field strengths (see, e.g. ref.  for a discussion of this point); in our simulations we will use

$`Z_{\mu \nu }`$ $`=`$ $`\mathrm{cos}\theta _\mathrm{W}n^a(x)W_{\mu \nu }^a\mathrm{sin}\theta _\mathrm{W}Y_{\mu \nu };`$
$`A_{\mu \nu }`$ $`=`$ $`\mathrm{sin}\theta _\mathrm{W}n^a(x)W_{\mu \nu }^a+\mathrm{cos}\theta _\mathrm{W}Y_{\mu \nu }.`$ (9)
Note that eqs. (7) to (9) reduce to the usual definitions away from the defect cores. We work in flat space and in the temporal gauge ($`W_0^a=Y_0=0`$ for $`a=1,2,3`$) so $`D_0\mathrm{\Phi }=_0\mathrm{\Phi }`$. With this gauge choice, Gauss’s Law becomes

$`_j(_0Y_j)`$ $`=`$ $`{\displaystyle \frac{i}{2}}\mathrm{sin}^2\theta _\mathrm{W}\left[\mathrm{\Phi }^{}_0\mathrm{\Phi }(_0\mathrm{\Phi })^{}\mathrm{\Phi }\right];`$
$`_j(_0W_j^a)ϵ^{abc}W_j^b_0W_j^c`$ $`=`$ $`{\displaystyle \frac{i}{2}}\mathrm{cos}^2\theta _\mathrm{W}\left[\mathrm{\Phi }^{}\tau ^a_0\mathrm{\Phi }(_0\mathrm{\Phi })^{}\tau ^a\mathrm{\Phi }\right],`$ (10)
(with $`j=1,2,3`$) which is then used to test the stability of the code.

As noted above, the $`\mathrm{sin}^2\theta _\mathrm{W}=1`$ case reduces to the semilocal case, where absolutely stable defects are known to exist. Setting one of the Higgs fields to zero further reduces this to the abelian Higgs model. We are able to use our knowledge of these systems, as well as the behaviour of artificially constructed infinite axially-symmetric $`Z`$-strings at various points in parameter space, to check the validity of our simulation code.

## 3 Numerical simulations

The equations of motion eq. (5) are discretized using a naïve staggered leapfrog method, where both the scalar and gauge fields are associated with lattice points. This procedure is the easiest to relate to previous work on semilocal defects , and we expect to have a fairly good picture at the scales we are interested in. A comparison to a link-variable discretization  of eq. (5) in the semilocal case ($`\theta _\mathrm{W}=\pi /2`$) shows that the changes from using a lattice implementation instead of this naïve one are well within the other uncertainties; see the appendix for a full analysis. The simulations are performed on a periodic cubic lattice whose time step is $`0.2`$ times the spatial step $`\mathrm{\Delta }t=0.2\mathrm{\Delta }x`$ ($`c=1`$), with *ad hoc* numerical viscosity terms added to each equation ($`\gamma \dot{\mathrm{\Phi }}`$, $`\gamma \dot{Y}`$ and $`\gamma \dot{W}^a`$ respectively) to reduce the system’s relaxation time. The expansion rate in an expanding universe would play this role, albeit as a time-dependent factor $`\gamma (t)`$ typically scaling as $`1/t`$. Several different values of $`\gamma `$ were tested, and generated very similar behaviour; throughout this work we use $`\gamma =0.5`$.

Two different strategies for setting the initial configurations were considered:

* Set all initial field velocities to zero and throw down random scalar field phases at every lattice point; average the field at each point with its 6 nearest neighbours and normalize the fields to the vacuum value ($`|\mathrm{\Phi }|^2=1`$ as $`\eta `$ has been rescaled out) iteratively (50 times), to get a smoother configuration; using these scalar field values, choose the initial values of the gauge fields to be

$`Y_\mu `$ $`=`$ $`0`$
$`W_\mu ^1`$ $`=`$ $`2\left(\psi _1_j\psi _4\psi _4_j\psi _1+\psi _3_j\psi _2\psi _2_j\psi _3\right)`$
$`W_\mu ^2`$ $`=`$ $`2\left(\psi _3_j\psi _1\psi _1_j\psi _3+\psi _4_j\psi _2\psi _2_j\psi _4\right)`$
$`W_\mu ^3`$ $`=`$ $`2\left(\psi _1_j\psi _2\psi _2_j\psi _1+\psi _4_j\psi _3\psi _3_j\psi _4\right)`$ (11)
where $`\mathrm{\Phi }^T=(\psi _1+i\psi _2,\psi _3+i\psi _4)`$, pseudo-minimizing the energy as in ref. . Note, however, that in the present case, the system does have enough gauge fields to cancel completelly the gradient energy. In fact, the configuration (11) sets to zero the gradient energy, the potential energy and the field-strength energy corresponding to the $`U(1)`$ gauge field, leaving the $`SU(2)`$ field-strength energy as the only non-zero term.
* Set all fields to zero, and also gauge field velocities to zero ($`\mathrm{\Phi }=W_\mu ^a=\dot{W}_\mu ^a=Y_\mu =\dot{Y}_\mu =0`$); provide some initial (smoothed) random velocities to the scalar field ($`\dot{\mathrm{\Phi }}0`$) following the general procedure as above. This initial configuration would be closer to that appropriate to a bubble nucleation scheme.

The overall results in our simulations did not show qualitatively different behaviour between the two cases. Indeed, in the first case the gauge field energy pseudo-minimization was ineffective enough that the scalar field typically began by climbing up the potential, restoring the symmetric phase, and rolled down again later. Previous experience with semilocal strings also shows that the results are fairly insensitive to the specific way initial conditions are implemented. For all these reasons, we chose the second initial configuration ($`b`$) for our investigations.

Interpreting simulations of electroweak string networks is more complicated than in the $`U(1)`$ cosmic string case. As in the semilocal case, electroweak strings are non-topological, and in this case there is no well-defined winding number, which makes identification of strings a more difficult task. We follow the strategy proposed in ref.  in order to study string formation: the system is evolved forward in time and we compute a set of gauge-invariant quantities at each time step, namely the $`Z`$\- and the $`A`$-field strengths given by eq. (9) and the modulus of the scalar field ($`|\mathrm{\Phi }|`$). String formation is then observed by visualizing isosurfaces in the $`Z`$-field strength ($`\sqrt{\frac{1}{2}Z_{ij}Z^{ij}}`$) and in the scalar field modulus.

After the rescaling eq. (4), it becomes clear that the only free parameters in our model are $`\beta `$ ($`\beta m_H^2/m_\mathrm{z}^2=8\lambda /g_\mathrm{z}^2`$) and $`\theta _\mathrm{W}`$. It is well known that with $`\theta _\mathrm{W}=\pi /2`$ and $`\beta <1`$ we have semilocal strings that are stable , and numerical simulations show segment formation and linkage . Moreover, for $`\mathrm{sin}^2\theta _\mathrm{W}1`$, $`\beta <1`$, there is a regime, albeit a rather narrow one, where infinitely long strings are perturbatively stable . Bearing these results in mind, the parameter space investigated in our simulations was $`0.9\mathrm{sin}^2\theta _\mathrm{W}1`$ and $`0.05\beta 1.5`$.

Simulation and visualization were performed on $`64^3`$, $`128^3`$ and $`256^3`$ lattices using the Cray T3E at NERSC, and high-performance computing facilities at University of Sussex and University of Groningen. Our quoted results all come from $`256^3`$ simulations.

## 4 Results

To test our code, we began by reproducing known results for cosmic and semilocal string networks, and for infinitely long, axially-symmetric $`Z`$-strings (made possible by our periodic boundary conditions as long as strings are simulated in pairs to keep the net flux equal to zero) in both the stable and unstable regimes, in particular verifying that the $`Z`$-strings disappeared in the unstable regime.

Having checked the code, we then ran it for values of the two free parameters ($`\theta _W,\beta `$) throughout the regime of interest, using the same initial conditions for each parameter pair. We chose to carry out one large simulation at each point in parameter space rather than many smaller ones, which gives improved dynamical range. As expected, after an initial transient, the typical configuration observed is a dumbbell — a segment of electroweak string joining a monopole/antimonopole pair. By continuity, for parameter values sufficiently close to the stable semilocal case we expect the monopole interactions to lead to the joining of strings to form longer segments.

The simulations show that some short segments of $`Z`$-string do join as expected. The joining rate is, however, lower than in the semilocal case, and decreases both as $`\mathrm{sin}^2\theta _\mathrm{W}`$ is decreased and as $`\beta `$ is increased. Thus, some segments which would eventually join in the semilocal case are seen instead to collapse in the electroweak case due to string tension. This is not surprising: in the semilocal case, we have global monopoles at the string ends, which have divergent scalar gradient energy and are more efficient at finding neighbouring monopoles. In the electroweak case, the monopoles at the string ends are proper magnetic monopoles, and the scalar gradients are cancelled much more efficiently by gauge fields. As $`\mathrm{sin}^2\theta _\mathrm{W}1`$ the cores become larger, and eventually overlap, making segments join. But as we move away from the semilocal case, the cores become smaller and the joining becomes less important.

Figure 1 shows timeslices of two typical simulations, with the different colours corresponding to the $`A`$\- and $`Z`$-field strengths.111Further colour images and movies can be found at
http://www.nersc.gov/$$borrill/defects/electroweak.html

The $`Z`$-field has a string-like form, whereas the $`A`$-field at the string ends is a spherical shell, corresponding to spherical magnetic monopoles. The $`A`$-field morphology can also be tube-like, denoting interaction between monopoles, illustrating the complexity of the overall dynamics. The interplay between string tension and monopole-antimonopole attraction causes some strings to shrink until they disappear, and others join to form longer strings.

To compare the number of defects in the semilocal and electroweak cases, we plot the number of lattice sites with $`Z`$-magnetic field strength ($`\sqrt{\frac{1}{2}Z_{ij}Z^{ij}}`$) greater than 25% of the maximum field strength found in the core of a NO string with the same value of $`\beta `$ (see figure 2). The number of lattice sites, and hence length of string, decreases more rapidly in the electroweak case as either $`\mathrm{sin}^2\theta _\mathrm{W}`$ decreases or $`\beta `$ grows.

At early times the configurations in different simulations are very similar. During the symmetry-breaking transient (the very first time steps) no defects exist. As the scalar field takes on a non-zero value a large number of very small string segments emerge, then the magnetic monopoles become visible at the ends of the segments due to the gathering of $`A`$-magnetic flux there. The upper panels of figure 1 show the $`Z`$\- and $`A`$-magnetic field strength of two different simulations at time $`t=50`$. Similarities in these first time steps, together with similarities using different initial conditions (see above), show that the initial configurations (i.e. the way the symmetry breaking is implemented) are not as important as the subsequent interaction between the scalar and gauge fields.

As the system evolves, we see that in only one of the two simulations do the small segments grow and connect to their neighbours, allowing the defect network to persist. In the lower panels the configuration at time $`t=200`$ can be seen. In one case there are long strings and connections are still happening, whereas in the other almost all of the dumbbells have annihilated. In particular, in the first case the final configuration contains strings which are much longer at the end of the simulation than any present in the early stages.

We wish to determine the parameter space for which the string network persists. The stability of the corresponding infinite $`Z`$-strings is a necessary condition for persistence, but not sufficient since the physically realizable collection of finite initial string segments may not link up. We therefore expect the persistence region of parameter space to lie entirely within the stability region.

To estimate the persistence of $`Z`$-strings we need to define a criterion to decide when defects are lasting long enough. There is clearly some degree of arbitrariness in how this criterion is set, and we have been guided by inspecting visually the evolution in different parameter regimes. In order to construct the persistence region in parameter space we simulated the system in $`256^3`$ boxes with periodic boundary conditions, and considered persistence to have occurred if at time $`t=200`$ there are more than 1000 lattice sites with a $`Z`$-magnetic flux greater than $`25\%`$ of the maximum of the NO simulation for the same $`\beta `$. For example, for $`\beta =0.3`$, we can see in figure 2 that persistence defined this way is exhibited only for $`\mathrm{sin}^2\theta _\mathrm{W}>0.995`$. From a suite of simulations using this criterion we obtain the ‘persistence limit’ shown in figure 3. As anticipated, the persistence region covers only a subset of the stability region, with persistence only for values of $`\mathrm{sin}^2\theta _\mathrm{W}`$ extremely close to one. Another possible criterion for persistence which can be easily automated is the presence at late times of strings that are reasonably long compared to their width.222To automate the calculation of the length of the individual strings, we compute the volume of each segment of string by counting the number of connected points whose $`Z`$-magnetic field is at least 25% of the maximum in a NO string at the same value of $`\beta `$. Then, the volume is divided by the cross-sectional area of the NO string. Finally, to obtain the length-to-width ratio, we divide it by the diameter of the corresponding NO string. The radius of the NO string is always taken to be the distance between the center and the point where the field strength is 25% of the maximum. At time $`t=200`$, the condition that there be strings at least five times their width leads to the second persistence line drawn in figure 3.

## 5 Conclusions

The numerical simulations described here show that a significant non-topological defect network can form in a generalized GSW model for an electroweak phase transition. The dynamics of such a network are extremely complicated, driven by string segment linkings and by isolated strings shrinking, and the details are highly sensitive to the two model parameters $`\theta _\mathrm{W}`$ and $`\beta `$. Our principal result is that in some regions of parameter space a persistent network of genuine non-topological defects can form. Though the actual version of electroweak theory in our Universe lies outside this parameter regime, these results show that, in models where topological (or semilocal) defects are possible in some limit, it is possible to get a network of non-topological defects close to that limit. Previous works in the literature  show that for parameters close to those permitting topological (or semilocal) defects, there is a regime where non-topological defects are stable. Our work confirms this result, and shows that, although in a narrower region, a sufficiently persistent network of defects can form in such a phase transition.

Note that the generation of the string segments is intrinsically dynamical and cannot be studied using initial condition arguments; in particular our results are compatible with those of . Here we are considering the time evolution of a network and not just the initial configuration. In fact, the first timesteps in such simulations correspond more to a numerical transient, in which physically reasonable initial conditions are established, than an actual phase transition. Only after this initial transient can the evolution of the network be trusted, and this is what determines whether the defects persist or decay. Given the primary role played by the gauge fields, it would be interesting to know if our conclusions generalize to non-topological string defects in other models such as the two-Higgs standard model .

###### Acknowledgments.

We thank Mark Hindmarsh and Tanmay Vachaspati for useful discussions, and the referees for their comments. AA and JU acknowledge support from grants CICYT AEN99-0315 and UPV 063.310-EB187/98. JU is also partially supported by a Marie Curie Fellowship of the European Community programme HUMAN POTENTIAL under contract number HPMT-CT-2000-00096. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098. JU is grateful to the Lawrence Berkeley National Laboratory at the University of California, the Kapteyn Institute and the Institute for Theoretical Physics at the University of Groningen, and the Astronomy Centre at the University of Sussex for their hospitality during visits, and the use of their computer facilities including those of the Sussex High-Performance Computing Initiative.

## Appendix A Discretization methods

In this work, we discretized the equations of motion eq. (5) replacing all the scalar and gauge fields by their values at the lattice points. For instance, a $`U(1)`$ covariant derivative was substituted by

$$D_j\varphi (x)=\frac{\varphi (x+x_j)\varphi (xx_j)}{2L}\frac{ig_\mathrm{Y}}{2}Y_i(x)\varphi (x)$$
(12)
where $`L`$ is the lattice spacing.

It is well known that this discretization scheme does not respect gauge invariance. Nevertheless, we are dealing with well-defined classical equations of motion in a particular gauge, and the breaking of gauge invariance should be irrelevant in this context. Indeed, we make use of the fact that Gauss’s Law is not automatically satisfied, and by monitoring it during the evolution we can check that our discretization is reasonably accurate. This discretization was chosen to facilitate comparison with earlier work .

There is an alternative discretization method widely used in the literature which protects gauge invariance and recovers the original equations of motion in the limit where the lattice spacing goes to zero . This method uses lattice link variables, i.e., the replacement of gauge fields by matrices living on the links between lattice points. In this case, a $`U(1)`$ covariant derivative will be substituted by

$$D_j\varphi (x)=\frac{1}{L}\left(e^{iLY_j(x)}\varphi (x+x_j)\varphi (x)\right).$$
(13)
In order to check whether our conclusions still hold using link variables, we first performed a series of simulations in the semilocal case to determine whether the evolution was compatible. Beginning from the same initial conditions, we evolved the system according to each discretization scheme. We observed that the simulations undergo extremely similar evolution on a pointwise basis, as seen in figure 4, though there are modest differences with slightly longer strings on the lattice gauge calculations, occasionally leading to extra connections.

To quantify this small difference and its impact on our statistical result, we performed a further 120 simulations using both discretizations. To save computing power we simulated our system in the semilocal case ($`\mathrm{sin}^2\theta _\mathrm{W}=1`$), discretized in $`64^3`$ cubes for different values of $`\beta `$. We used several different initial conditions, but for each one we let the system evolve using both the naïve discretization of the (lagrangian) equations of motion and the link variable method for discretizing the (hamiltonian) equations of motion. The initial conditions were calculated using the method b) described in the text (section 3), and as in the electroweak simulations performed in this work, we added an *ad hoc* damping term ($`\eta =0.5`$). As our interest lies mainly in the late-time behaviour of the dumbbell network, because the smaller cube size forces $`t64`$, we chose to calculate the total string length at $`t=56`$ for both cases.

Figure 5 shows the results. We computed the string lengths as explained in section 4, footnote 2, and then added all string lengths for strings longer than five times their length. That final measure is represented in figure 5, which shows clearly that the differences between the schemes are well within the uncertainties. Although performed for only $`64^3`$ cubes, this result should continue to hold for the larger cubes used to derive the main results in the body of this paper. The system simulated in this test is only the semilocal case, not the full electroweak case, but we are convinced that using the lattice link variable method will not alter significantly the results presented in this paper.


text:
# All Optical Formation of an Atomic Bose-Einstein Condensate

## Abstract

We have created a Bose-Einstein condensate of $`{}_{}{}^{87}\text{Rb}`$ atoms directly in an optical trap. We employ a quasi-electrostatic dipole force trap formed by two crossed $`\text{CO}_2`$ laser beams. Loading directly from a sub-doppler laser-cooled cloud of atoms results in initial phase space densities of $`1/200`$. Evaporatively cooling through the BEC transition is achieved by lowering the power in the trapping beams over $`2`$ s. The resulting condensates are $`F=1`$ spinors with $`3.5\text{ x }10^4`$ atoms distributed between the $`m_F=(1,0,1)`$ states.

03.75.Fi, 32.80.Pj, 03.67.-a

The first observation of Bose Einstein condensates (BEC) in dilute atomic vapors in a remarkable series of experiments in 1995  has stimulated a tremendous volume of experimental and theoretical work in this field. Condensates are now routinely created in over 30 laboratories around the world, and the pace of theoretical progress is equally impressive . The recipe for forming a BEC is by now well-established . The atomic vapor is first pre-cooled, typically by laser cooling techniques, to sub-mK temperatures and then transferred to a magnetic trap. Further cooling to BEC is then achieved by evaporatively cooling the atoms in the magnetic trap using energetically selective spin transitions .

All-optical methods of reaching the BEC phase transition have been pursued since the early days of laser cooling. Despite many impressive developments beyond the limits set by Doppler cooling, including polarization gradient cooling , VSCPT , Raman cooling , and evaporative cooling in optical dipole force traps , the best efforts to-date produce atomic phase space densities $`n\lambda _{dB}^3`$ a factor of 30 away from the BEC transition . The principal roadblocks have been attributed to density-dependent heating and losses in laser cooling techniques, residual heating in optical dipole force traps or the unfavorable starting conditions for evaporative cooling. Hence, optical traps have played only an ancillary role in BEC experiments. The MIT group used a magnetic trap with an ‘optical dimple’ to reversibly condense a magnetically confined cloud of atoms evaporatively cooled to just above the phase transition . Additionally, Bose condensates created in magnetic traps have been successfully transferred to shallow optical traps for further study . In all these cases, however, magnetic traps provided the principle increase of phase space density (by factors up to $`10^6`$) to the BEC transition.

In this letter, we present an experiment in which we have created a Bose condensate of 87Rb atoms directly in an optical trap formed by tightly focused laser beams. Following initial loading from a laser cooled gas, evaporative cooling through the BEC transition is achieved by simply lowering the depth of the optical trap. Our success is due in part to a high initial phase space density realized in the loading of our optical dipole trap and in part to the tight confinement of the atoms that permits rapid evaporative cooling to the BEC transition in $`2`$ s. This fast evaporation relaxes considerably the requirement for extremely long trap lifetimes typical of magnetic trap BEC experiments. Additionally, in contrast to magnetic traps which only confine one or two magnetic spin projections , our technique is spin-independent and the condensates that we form are $`F=1`$ three-component spinors .

We utilize a crossed-beam optical dipole force trap employing tightly focused high-powered (12 W) infrared ($`\lambda =10.6`$ $`\mu `$m) laser beams. For trapping fields at this wavelength, the trapping potential for the ground state atoms is very well approximated by $`U(𝐫)=\frac{1}{2}\alpha _g\left|E\left(𝐫\right)\right|^2`$ where $`\alpha _g`$ is the ground state dc-polarizability of the atom (5.3 x 10-39 m2C$`/`$V for atomic rubidium) and $`E\left(𝐫\right)`$ is the electric field amplitude. A significant feature of these ‘quasielectrostatic’ traps (QUESTs)  is that heating due to spontaneous emission of the atoms is completely negligible. In contrast to previous QUESTs that have employed single focused beams  or standing wave configurations , we employ here a cross-beam geometry  to provide a balance of tight confinement in three dimensions ($`1.5`$ kHz oscillation frequencies at full power) and a relatively large loading volume.

Our experiments begin with a standard vapor loaded MOT. The trapping beams consist of three orthogonal retroreflected beams in the $`\sigma ^+\sigma ^{}`$ configuration. They are tuned 15 MHz below the 5S1/2-5P3/2 $`F=2F^{}=3`$ transition of 87Rb and each of the three beams has a waist of 0.7 cm and a power of 25 mW. A repump beam tuned to the $`F=1F^{}=2`$ transition is overlapped with one of the MOT trapping beams. The MOT is loaded for 5 s directly from the thermal vapor during which we collect 30 x 106 atoms. After loading the MOT, the trap configuration is changed to maximize the transfer of atoms to the optical trap. The repump intensity is first lowered to $``$10 $`\mu `$W/cm2 for a duration of 20 ms, and then the MOT trap beams are shifted to the red by 140 MHz for a duration of 40 ms. At this point the MOT beams are extinguished and the current in the MOT coils is turned off. In order to optically pump the atoms into the $`F=1`$ hyperfine states, the repump light is shuttered off 1 ms before the trap beams are extinguished; we measure the efficiency of the optical pumping to the $`F=1`$ state to be $`>95\%`$. The CO2 laser beams are left on at full power (12 W) throughout the MOT loading and dipole trap loading process.

The trapping beams are generated from a commercial CO2 gas laser ($`\lambda =10.6`$ $`\mu `$m). The beams are tightly focused and intersected at right angles; one beam is oriented in the horizontal direction and one beam is inclined at 45$`^\text{o}`$ from the vertical direction. Each beam passes through an acousto-optic modulator to provide independent control of the power in the two beams. Additionally, the beams are frequency shifted 80 MHz relative to each other so that any spatial interference patterns between the two beams are time-averaged to zero . Each beam has a maximum power of 12 W, and the beams are focused to a minimum waist $``$ 50 $`\mu `$m with $`f=38`$ mm focal length ZeZn aspherical lenses inside the chamber.

Following standard techniques, the number of trapped atoms and their momentum distribution are observed using absorptive imaging of the released atoms. The trapped atoms are released by suddenly ($`<`$ 1 $`\mu `$s) switching off the trapping laser beams. Following a variable ballistic expansion time (typically 2-20 ms), the cloud is illuminated with a 50 $`\mu `$s pulse of $`F=1F^{}=2`$ light applied to the atoms concurrent with a vertically oriented circularly polarized probe beam tuned to the $`F=2F^{}=3`$ transition. The probe intensity is $``$0.3 mW/cm2, and each atom scatters up to 150 photons from the probe with no observable blurring of the cloud. The shadow of the atom cloud is imaged onto a slow-scan CCD camera. The measured spatial resolution of our imaging system is $`10`$ $`\mu `$m. From these images, the number of atoms and their temperature are determined. Together with the trap oscillation frequencies (which are measured directly using parametric excitation) the spatial density, elastic collision rate, and phase space density can be derived.

We first discuss the properties of the trap following loading from the MOT without employing forced evaporative cooling. The earliest that we observe the trapped atoms is 100 ms after loading to allow for the untrapped atoms to fall away. At this time, the trap contains 2 x 106 atoms at a temperature of 75 $`\mu `$K. Maintaining full power in the trap beams, we observe a rapid evaporation from the trap in the first 1.5 s, during which of two-thirds of the trapped atoms are lost and the temperature falls to 38 $`\mu `$K. The relative phase space density increases by a factor of 3 during this stage. This rapid evaporation gives way to a much slower exponential trap decay with a time constant of 6(1) s. The temperature continues to fall gradually to a final value of 22 $`\mu `$K on a 10 s timescale.

The mean frequency of the trap at these powers is measured to be 1.5(2) kHz, from which we infer an initial peak density of the trap of $``$2 x 1014 atoms$`/`$cm$`^3.`$ The initial phase space density following the rapid free evaporation stage is calculated to be 1/200 (or 1/600 assuming equal distribution of the atoms in the three trapped internal states). Using the measured value for the three body loss rate  yields an initial 3-body loss rate constant of 1.3 s-1. Hence it is possible that 3-body inelastic collisions contribute to the rapid initial loss of atoms, however we note that such collisions do not explain the observed cooling. The estimated elastic collision rate at this point is a remarkable 12 x 103 s-1—higher than the trap oscillation frequency. Although the derived quantities must be taken with some caution, we note that densities of 3 x 1013 atoms$`/`$cm-3 and a phase space density of 1/300 have been reported in a 1-D CO2 lattice trap with 85Rb (albeit with much fewer atoms) , and we have measured similar results in our laboratory in 1-D CO2 lattice traps with 87Rb. It is not clear what yields these extremely high initial densities in these traps at these low temperatures. A blue detuned Sisyphus process is suggested in  (densities of $``$1013 atoms$`/`$cm-3 have been observed in a dipole trap with blue detuned Sisyphus cooling ), but we suspect that there is an effective dark SPOT  in the repump beam at the trap center induced by the AC stark shift due to the dipole trap beams, which reduces density-limiting interactions with the MOT light fields. Using our trap parameters, we estimate that the repump scattering rate is reduced at the trap center by a factor of 40 relative to outside the trap. Our improved performance even relative to the 1-D lattice traps may be due to more efficient loading from the ‘tails’ of the crossed trap geometry. In any case, the net result for our trap is a very high initial phase space density combined with a large number of atoms and a fast thermalization time—together these conditions provide a very favorable starting point for evaporative cooling.

Efficient forced evaporation requires selectively ejecting the more energetic atoms from the trap such the the remaining atoms rethermalize at a lower temperature with a higher net phase space density. For optical traps, the simplest way to force evaporation is to lower the trap depth by decreasing the power in the trap beams. This technique was used in one of the first demonstrations of evaporative cooling in alkali atoms , where, starting with only 5000 atoms, a phase space density increase of $`30`$ was realized. The drawback to this method is that lowering the trap depth also lowers the trap oscillation frequency and the rethermalization rate. Hence the evaporation rate can slow down prohibitively. In our case, although the rethermalization rate falls by a factor of 50 by the end of the evaporation cycle, it nonetheless remains fast enough to allow us to reach the BEC transition in 2.5 s as described below.

Our procedure for forced evaporation is as follows: immediately after loading the trap, the power in both trap beams is ramped to 1 W in 1 s. The power in both beams is then ramped to a variable final power in 1 s and maintained at this low level for 0.5 s, after which the atoms are released from the trap and imaged as discussed above. Fig. 1 shows three such images at different final powers for the trap lasers. The left image shows a cloud well above the BEC transition point, and reveals an isotropic gaussian momentum distribution expected for a thermal cloud of atoms in equilibrium. As the evaporative cooling proceeds to lower powers, a bimodal momentum distribution appears with a central non-isotropic component characteristic of Bose condensates (Fig. 1, center). As the power is lowered further, the central, non-spherical peak of the distribution becomes more prominent, and the spherical pedestal diminishes. At a trap power of 190 mW (Fig. 1, right), the resulting cloud is almost a pure condensate and contains 3.5 x 104 atoms. For lower powers, the cloud rapidly diminishes as the trap can no longer support gravity. Line profiles of the images in Fig. 1 are shown in Fig. 2 for quantitative comparison. The profiles are taken along the orientation of the minor axis of the condensate. Also shown are gaussian fits to the data in the wings that clearly show the bimodal nature of the momentum distributions near the BEC phase transition.

The critical condensation temperature $`T_c`$ is a function of the number of atoms $`N`$ and the mean trap frequencies $`\overline{\omega }=(\omega _1\omega _2\omega _3)^{1/3}`$ according to $`k_BT_c=\mathrm{}\overline{\omega }(N/1.202)^{1/3}.`$ Below the critical temperature, the condensate fraction should grow as $`N_0/N=1(T/T_c)^3`$. We obtain the condensate fraction and the temperature of the normal component from 2D fits to the absorptive images using the methods described in . In Fig. 3, the condensate fraction is plotted versus normalized temperature, $`T/T_c(N,\overline{\omega })`$ near the critical point, where we assume that $`\overline{\omega }^2`$ is proportional to the trap power. The agreement with the theoretical curve is reasonable given the scatter in the data. At the trap power of 350 mW, which is near the critical point, the trap contains 180,000 atoms at a temperature of 375 nK. Together with the measured trap frequencies at this power $`\omega _1,\omega _2,\omega _3=2\pi (72,175,350)`$ Hz, we infer a phase space density at this point of 1.4 (or 0.45 if the three internal states are equally populated), a density of 4.8 x $`10^{13}`$ cm-3 and an elastic collision rate of 300 s-1.

We have measured the growth of the freely expanding condensate for expansion times from 5-20 ms. The measured aspect ratio of the cloud is 1.5(0.1) and shows a slight increase for longer expansion times. The measured size of the minor axis of the cloud grows linearly with time (within the limits of our imaging resolution) to a size of 100 $`\mu `$m (full-width) at 15 ms. The observed 2D projection of the cloud is oriented at a $`30^\text{o}`$ angle relative to the horizontal trap. Quantitative comparison of the observed condensates with theory requires knowledge of the trap oscillation frequencies and as well as orientation of the principal axes of the trap. Although we can measure the frequencies, the principal axes depend on the beam ellipticity and alignment. In our case, the trap beams are aberrated due to off-center propagation through the focussing lenses (required to overlap the beams), and the fact that the laser output beam is not a pure gaussian TEM00 spatial mode. Nonetheless, using the measured trap frequencies, scaled for power, mean field theory (see e.g. ) predicts aspect ratios of 1.2-1.7 for a 15 ms expansion and cloud sizes of 60-130 $`\mu `$m, consistent with our observations.

The $`1/e`$ lifetime of the condensate is measured to be 3.5(1) s, which is somewhat smaller than the 6(1) s lifetime observed in the trap. We observe no residual heating of the condensate for the lifetime of the condensate, although because the trap potential is quite shallow, it is possible that any heating would be offset by subsequent evaporation from the trap.

The trapped atoms are optically pumped during loading into the $`F=1`$ groundstate, however, no attempt is made to further optically pump the atoms into a single $`m_F`$ state. Hence we expect the population to be a mixture of the $`m_F=1,0,1`$ spin projections. To measure the spin content of our condensate, we apply a weak field gradient to the cloud after the cloud is released from the dipole trap . An absorptive image of the cloud is shown in Fig. 4, and reveals that the condensate is comprised of a mixture of spin-states. Each cloud has a similar non-spherical momentum distribution characteristic of condensates. Determining the relative weights of the spin states is complicated by the differing matrix elements of the probe transitions. Nonetheless, if we assume that the atoms are all optically pumped by the circularly polarized probe and that all atoms scatter the same number of photons, we infer a 3:1:1 weighting of the $`m_F=1,0,1`$ mixture respectively for the image shown in Fig. 4.

In summary, we have realized a Bose condensate of $`{}_{}{}^{87}\text{Rb}`$ atoms directly in an optical trap. This technique seems promising for creating and studying more complex spinor condensates (e.g. $`{}_{}{}^{85}\text{Rb},F=2`$) as well as multi-atom mixtures. Additionally, by eliminating the need for a magnetic trap, it may be possible to realize condensates in atoms or molecules lacking a suitable magnetic moment, perhaps by using sympathetic cooling to precool the sample. Finally, our technique offers considerable experimental simplicity and speed, easing the requirement for ultra high vacuum environments, and eliminating the need for strong magnetic trapping fields.

We would like to acknowledge the technical assistance of D. Zhu and helpful discussions with T.A.B. Kennedy, C. Raman, S. Yi, and L. You. This work was supported by the National Security Agency (NSA) and Advanced Research and Development Activity (ARDA) under Army Research Office (ARO) contract number DAA55-98-1-0370.


text:
# Einstein gravity as an emergent phenomenon?

## 1 Analog gravity

There is a risk that all current attempts at “quantizing gravity” are condemned to failure because they have been starting from fundamentally flawed premise, and that in reality there are no fundamental gravitational degrees of freedom to quantize — it is possible that Einstein gravity is an emergent phenomenon, in the same sense that fluid dynamics emerges from molecular physics as a low-momentum long-distance approximation. In this essay we will take a careful look at this idea, and highlight some of the possibilities, problems, and opportunities that such a situation entails.

We were led to these notions via current research on analog models of general relativity . Because of the extreme difficulty (and inadvisability) of working with intense gravitational fields in a laboratory setting, interest has now turned to investigating the possibility of *simulating* aspects of general relativity — though it is not *a priori* expected that all features of Einstein gravity can successfully be carried over to the analog models. Numerous rather different physical systems have now been seen to be useful for developing analog models of general relativity. A literature search as of March 2001 finds well over a hundred scientific articles devoted to one or another aspect of analog gravity and effective metric techniques. The sheer number of different physical situations lending themselves to an “effective metric” description strongly suggests that there is something deep and fundamental going on.

Typically these are models of general relativity, in the sense that they provide an effective metric and so generate the basic kinematical background in which general relativity resides; in the absence of any dynamics for that effective metric we cannot really speak about these systems as models for general relativity. However, as we will discuss more fully bellow, quantum effects in these analog models might provide of some sort of dynamics resembling general relativity.

Remember that for mechanical systems with a finite number of degrees of freedom small oscillations can always be resolved into normal modes: a finite collection of uncoupled harmonic oscillators. For a classical field theory you would also expect similar behaviour: small deviations from a background solution of the field equations will be resolved into travelling waves; then these travelling waves can be viewed as an infinite collection of harmonic oscillators, or a finite number if the field theory is truncated in the infra-red and ultra-violet, to which you can then apply a normal mode analysis. The physically interesting question is whether this normal mode analysis for field theories can then be reinterpreted in a “geometrically clean” way in terms of some “effective metric” and “effective geometry”. In many cases the answer is definitely yes: Linearization of a Lagrangian-based dynamics, or linearization of any hyperbolic second-order PDE, will automatically lead to an effective Lorentzian geometry that governs the propagation of the fluctuations. \[The most general situation (multiple scalar fields, or a multi-component vector or tensor) is quite algebraically messy — details of that situation will be deferred for now.\]

Once the notion of a derived “effective metric” has been established, we can certainly consider the effect of quantizing the linearized fluctuations. At one loop the quantum effective action will contain a term proportional to the Einstein–Hilbert action — this is a key portion of Sakharov’s “induced gravity” idea . In the closing segment of the essay we argue that the occurrence of not just an “effective metric”, but also an “effective geometrodynamics” closely related to Einstein gravity, is a largely unavoidable feature of the linearization and quantization process.

## 2 Effective metric

Suppose we have a single scalar field $`\varphi `$ whose dynamics is governed by some first-order Lagrangian $`(_\mu \varphi ,\varphi )`$. (By “first-order” we mean that the Lagrangian is some arbitrary function of the field and its first derivatives.) We want to consider linearized fluctuations around some background solution $`\varphi _0(t,\stackrel{}{x})`$ of the equations of motion, and to this end we write

$$\varphi (t,\stackrel{}{x})=\varphi _0(t,\stackrel{}{x})+ϵ\varphi _1(t,\stackrel{}{x})+O(ϵ^2).$$
(1)
Linearizing the Euler–Lagrange equations results in a second-order differential equation with position-dependent coefficients (these coefficients all being implicit functions of the background field $`\varphi _0`$). Following an analysis developed for acoustic geometries (Unruh , Visser *et al* ), which also applies to this much more general situation, this can be given a nice clean geometrical interpretation in terms of a d’Alembertian wave equation — provided we *define* the effective spacetime metric by

$$\sqrt{g}g^{\mu \nu }f^{\mu \nu }\left\{\frac{^2}{(_\mu \varphi )(_\nu \varphi )}\right\}|_{\varphi _0},$$
(2)
that is,

$$g_{\mu \nu }(\varphi _0)=\left(\mathrm{det}\left\{\frac{^2}{(_\mu \varphi )(_\nu \varphi )}\right\}\right)^{1/(d1)}|_{\varphi _0}\left\{\frac{^2}{(_\mu \varphi )(_\nu \varphi )}\right\}^1|_{\varphi _0}.$$
(3)
The equation of motion for the linearized fluctuations can then be written in the geometrical form

$$\left[\mathrm{\Delta }(g(\varphi _0))V(\varphi _0)\right]\varphi _1=0,$$
(4)
where $`\mathrm{\Delta }`$ is the d’Alembertian operator associated with the effective metric $`g(\varphi _0)`$, and $`V(\varphi _0)`$ is a background-field-dependent potential. It is important to realise just how general the result is: it works for any Lagrangian depending only on a single scalar field and its first derivatives. The linearized PDE will be *hyperbolic* (and so the linearized equations will have wave-like solutions) if and only if the effective metric $`g_{\mu \nu }`$ has Lorentzian signature $`\pm [,(+)^d]`$ . Note that $`d=1`$ space dimensions is special, and the present formulation does not work unless $`\mathrm{det}(f^{\mu \nu })=1`$. This observation can be traced back to the conformal covariance of the Laplacian in $`1+1`$ dimensions, and implies (perhaps ironically) that the only time the procedure risks failure is when considering a field theory defined on the world sheet of a string-like object.

Indeed, even if you do not have a Lagrangian, it is still possible to extract an “effective metric” for a system with one degree of freedom. (More precisely, we can define a conformal class of effective metrics. The analysis is not as geometrically “clean”.) While several of the technical details are different from the Lagrangian-based analysis, the basic flavor is the same: The key point is that hyperbolicity of the linearized PDE is defined in terms of the presence of a matrix of indefinite signature $`\pm [,(+)^d]`$ . This matrix is enough to define a conformal class of Lorentzian metrics, and picking the “right” member of the conformal class is largely a matter of taste — do whatever makes the “geometrized” equation look cleanest.

## 3 Effective dynamics

At this stage we have derived the existence of a background metric $`g_{\mu \nu }(\varphi _0)`$ and linearized fluctuations governed by the equation (4). We shall now try to see how and to what extent it is possible to define a *dynamics* for this metric. Of course it would be not particularly useful to search for such dynamics in the equations of motion of the “fundamental system” \[those derived from $`(\varphi ,\varphi )`$\]. It is from this perspective obvious that the dynamical equations of the effective metric should also be regarded as an “emergent” phenomenon. The idea of “induced gravity”, proposed several years ago by Andrei Sakharov , provides a natural framework for such an emergence of an effective geometrodynamics. Along these lines we shall now try to derive Einstein-like equations from the one-loop action of the $`\varphi `$ field.

Expand the field $`\varphi `$ into a background $`\varphi _b`$, which does not necessarily satisfy any classical equation of motion, and a quantum fluctuation $`\varphi _q`$ so that $`\varphi =\varphi _b+\varphi _q`$. Integrate out the quantum fluctuations; then at one loop

$$\mathrm{\Gamma }[g(\varphi _b),\varphi _b]=S[\varphi _b]+\frac{1}{2}\mathrm{}\mathrm{tr}\mathrm{ln}\left[\mathrm{\Delta }(g(\varphi _b))V(\varphi _b)\right]+O(\mathrm{}^2).$$
(5)
Here the determinant of the differential operator may be defined in terms of zeta functions or heat kernel expansions . Note also that the effective action depends on the background field in two ways: *explicitly* through $`\varphi _b`$, and *implicitly* through $`g(\varphi _b)`$. The key point is that defining the determinant requires both regularization and renormalization, and that doing so introduces counterterms proportional to the first $`d/2`$ Seeley-DeWitt coefficients . The zeroth Seeley–DeWitt coefficient $`a_0`$ induces a cosmological constant, while $`a_1`$ induces an Einstein–Hilbert term, and there are additional terms proportional to $`a_2`$. All in all:

$`\mathrm{\Gamma }[g(\varphi _b),\varphi _b]`$ $`=`$ $`S[\varphi _b]+\mathrm{}{\displaystyle \sqrt{g}\kappa \left[2\mathrm{\Lambda }+R(g)6V(\varphi _b)\right]\mathrm{d}^{d+1}x}`$ (6)
$`+\mathrm{}X[g(\varphi _b),\varphi _b]+O(\mathrm{}^2).`$
Here $`X[g(\varphi _b),\varphi _b]`$ denotes all other finite contributions to the renormalized one-loop effective action. It is the automatic emergence of the Einstein–Hilbert action as part of the one-loop effective action that is the salient point. Note that our approach is not identical to Sakharov’s idea — in his proposal the metric was put in by fiat, but without any intrinsic dynamics; all the dynamics was generated via one loop quantum effects. In our proposal the very existence of the effective metric itself is an emergent phenomenon. In Sakharov’s approach the metric was free to be varied at will, leading precisely to the Einstein equations (plus quantum corrections); in our approach the metric is not a free variable and the equations of motion will be a little trickier.

The quantum equations of motion are defined in the usual way by varying $`\mathrm{\Gamma }[\varphi _b]`$ with respect to the background $`\varphi _b`$. It is important to remember that the metric is a function of the background field so that it does not make sense to vary the metric independently — we must always evaluate variations using the chain rule.

$$\frac{\delta \mathrm{\Gamma }[g(\varphi _b),\varphi _b]}{\delta \varphi _b(x)}\frac{\delta \mathrm{\Gamma }[g(\varphi _b),\varphi _b]}{\delta \varphi _b(x)}|_{g_b}+\frac{\delta \mathrm{\Gamma }[g(\varphi _0),\varphi _0]}{\delta g_{\mu \nu }}|_{\varphi _b}\frac{\delta g_{\mu \nu }(\varphi _b)}{\delta \varphi _b(x)}.$$
(7)
In this way we see that the equations of motion have a part coming from the variation with respect to the background field plus a part proportional to the variation with respect to the metric. It is this second part which provides an Einstein-like dynamics. The presence of the terms generated by variation with respect to $`\varphi _b`$ leads to the interesting conclusion that one needs to assume, in order to get a dynamics as close as possible to that of Einstein, that these terms satisfy a special constraint: there should exist a functional $`Y[g]`$, depending only on the effective metric, such that

$`\left\{{\displaystyle \frac{\delta S[\varphi _b]}{\delta \varphi _b}}6\mathrm{}\kappa {\displaystyle \frac{\delta \sqrt{g}V(\varphi _b)}{\delta \varphi _b}}+\mathrm{}{\displaystyle \frac{\delta X[g(\varphi _b),\varphi _b]}{\delta \varphi _b}}|_{g_b}\right\}`$ $`=`$ $`\mathrm{}{\displaystyle \frac{\delta Y[g]}{\delta g_{\mu \nu }}}{\displaystyle \frac{\delta g_{\mu \nu }(\varphi _b)}{\delta \varphi _b(x)}}`$ (8)
$`+O(\mathrm{}^2).`$
In this case the background geometry decouples from the effective metric and we have

$$\left[\kappa \left(G^{\mu \nu }(g)+\mathrm{\Lambda }g^{\mu \nu }\right)+\frac{1}{\sqrt{g}}\frac{\delta \{X[g(\varphi _b)]+Y[g(\varphi _b)]\}}{\delta g_{\mu \nu }}\right]\frac{\delta g_{\mu \nu }(\varphi _b)}{\delta \varphi _b(x)}=O(\mathrm{}).$$
(9)
The $`\delta X[g(\varphi _b)]/\delta g`$ term denotes the type of “curvature squared” correction to the Einstein equations that is commonly encountered in string theory (indeed in almost any candidate theory for quantum gravity), and also in the usual implementation of Sakharov’s approach. However it must be emphasised that because of the contraction with the $`\delta g_{\mu \nu }(\varphi _b)/\delta \varphi _b(x)`$ these are not the usual Einstein equations, though they are certainly implied by the (curvature enhanced) Einstein equations. It is in this sense that we can begin to see the structure of Einstein gravity emerging from this field-theoretic normal mode analysis.

## 4 Prospects

In this essay we have argued that the emergence of an “effective metric”, in the sense that this notion is used in the so-called “analog models” of general relativity, is a rather generic feature of the linearization process. While the existence of an effective metric by itself does not allow you to simulate all of Einstein gravity, it allows one to do quite enough to be really significant — in particular it seems that the existence of an effective Lorentzian metric is really all that is in principle needed to obtain simulations of the Hawking radiation effect .

By invoking one-loop quantum effects, we have argued that something akin to Sakharov’s induced gravity scenario is operative: in particular we can generically argue that there is a term in the quantum effective action proportional to the Einstein–Hilbert action. However because of the technical assumption that the effective metric depends on the background only via the scalar field $`\varphi _0(x)`$ we have not been able to reproduce full Einstein gravity, though certainly have some extremely suggestive results along this line.

The major steps that are needed to extend this idea to a full fledged theory are

1. The question of what happens when many fields are present in the problem: the major piece of additional physics is the possible presence of birefringence, or more generally “multi-refringence”, with different normal modes possibly reacting to different metrics. The Eötvös experiment \[the observational universality of free fall to extremely high accuracy\] indicates that all the physical fields coupling to ordinary bulk matter “see” to high precision the *same* metric, allowing us to formulate the Einstein Equivalence principle and speak of *the* metric of spacetime.
2. Whether the addition of extra fields helps one to obtain a better approximation to full Einstein gravity — this because you would get one equation of motion per background field, so with six or more fields you would expect to be able to explore the full algebraic structure of the metric. So adding extra fields, which is technically a hindrance in the kinematical part of the program (developing the effective metric formalism), should in compensation allow one to more closely approach the dynamics of Einstein gravity.

In summary: The full generality of the situations under which effective metrics are encountered is truly remarkable, and the extent to which the resulting analog models seem able to reproduce key aspects of Einstein gravity is even more remarkable. The physics of these systems is fascinating, and the potential for laboratory investigation of models close to (but not necessarily identical to) Einstein gravity is extremely encouraging.

Our interpretation of these results is that they provide suggestive evidence that what we call Einstein gravity (general relativity) is an almost automatic low-energy consequence of almost any well behaved quantum field theory: the occurrence of an effective metric is almost automatic (even in the classical theory), while the presence of Einstein-like dynamics is almost guaranteed by one-loop quantum effects.

## Acknowledgments

The authors wish to thank Sebastiano Sonego for his comments and suggestions.