text:
# Vortex phase diagram in rotating two-component Bose-Einstein condensates

## Abstract

We investigate the structure of vortex states in rotating two-component Bose-Einstein condensates with equal intracomponent but varying intercomponent coupling constants. A phase diagram in the intercomponent-coupling versus rotation-frequency plane reveals rich equilibrium structures of vortex states. As the ratio of intercomponent to intracomponent couplings increases, the interlocked vortex lattices undergo phase transitions from triangular to square, to double-core lattices, and eventually develop interwoven “serpentine” vortex sheets with each component made up of chains of singly quantized vortices.

A system described by a multicomponent order parameters can excite various exotic topological defects. There have been great interest in the study of such topological defects, which has been explored extensivly in condensed matter systems such as anisotropic superfluid 3He helium3 . This study is also closely related with some problems in superconductors of UPt3 and Sr2RuO4 which are thought to have non-s-wave symmetry supercond  and high-energy physics and cosmology vorbook . Dilute atomic Bose-Einstein condensates (BECs) offer another superior testing ground to investigate the physics of exotic topological defects because they are free from impurities and well controlled by optical techniques. A trapped BEC in a rotating potential generates quantized vortices which are topological defects characteristic of superfluidity, and several groups have suceeded in observing ordered vortex lattices Madison . While many interesting phenomena have been found in single-component BECs Madison , a rich variety of static and dynamic phenomena are expected to occur in a system of rotating two-component BECs consisting, for example, of two different hyperfine spin states of atoms Hall . This system is characterized by three coupling constants denoted by $`C_{11}`$ and $`C_{22}`$ (for intracomponents) and $`C_{12}`$ (for intercomponent) Ho ; Ersy ; Ohberg1 ; Pu ; Tim . In this Letter, we report equilibrium structures of vortex states in two-component BECs with $`C_{11}=C_{22}C`$ but varying values of $`C_{12}`$ and the rotation frequency based on the numerical analysis of the Gross-Pitaevskii equation.

Mueller and Ho studied the vortex lattice structure of two-component BECs by assuming the lowest Landau level approximation and a perfect lattice Mueller . They proposed a phase diagram on the lattice structure and found that as $`C_{12}/C`$ is increased, the system undergoes a continuous structural change from triangular to square lattices. Kita et al. discussed the lattice structure in $`F=1`$ spinor BECs Kita . However, both studies assume the fast rotation limit where the rotation frequency is very close to the radial trapping frequency, so a vast experimentally accessible parameter regime remains to be explored. When there is no rotation, it is known that we obtain miscible condensates as the ground state for $`C_{12}/C<1`$, while for $`C_{12}/C>1`$ the two components phase separate Ho ; Ersy ; Tim . In the presence of rotation, our numerical analysis reveals a rich variety of vortex states which have eluded analytic treatments. In particular, we find that for $`C_{12}/C1`$ double-core vortex lattices or interwoven “serpentine” vortex sheets become stabilized for a considerable range of rotation frequency; these states can therefore be observed in current experimental situations.

The equilibrium solutions of two-component BECs in a frame rotating with the angular velocity $`𝛀=\mathrm{\Omega }\widehat{𝐳}`$ are obtained by solving the time-independent coupled Gross-Pitaevskii (GP) equations for the condensate wave functions $`(h_i+_{j=1,2}C_{ij}|\mathrm{\Psi }_j|^2)\mathrm{\Psi }_i=\mu _i\mathrm{\Psi }_i`$ ($`i=1,2`$) Ho ; Ersy ; Ohberg1 ; Pu ; Tim , where $`h_i=(m_{12}/m_i)^2+V_i(\mathrm{\Omega }/\overline{\omega })L_z`$ with the reduced mass $`m_{12}=m_1m_2/(m_1+m_2)`$ and the harmonic potential $`V_i=(m_i\omega _i^2/4m_{12}\overline{\omega }^2)r^2`$ . Here, length and energy are measured in units of $`b_{\mathrm{ho}}=\sqrt{\mathrm{}/2m_{12}\overline{\omega }}`$ and $`\mathrm{}\overline{\omega }\mathrm{}(\omega _1+\omega _2)/2`$, respectively. In this work, the numerical calculation is done in the two-dimensional $`xy`$ space. The wave functions are normalized as $`𝑑x𝑑y|\mathrm{\Psi }_i(x,y)|^2=1`$ and, then, the intra- and intercomponent coupling constants are written as $`C_{ii}=8\pi (m_{12}/m_i)N_i^{2\mathrm{d}}a_{ii}`$ and $`C_{ij}=4\pi N_j^{2\mathrm{d}}a_{12}`$ $`(ij)`$ with the corresponding s-wave scattering lengths $`a_{11}`$, $`a_{22}`$, $`a_{12}`$ (assumed to be positive), and the particle numbers $`N_i^{2\mathrm{d}}`$ per unit length along the $`z`$-axis.

For simplicity, the number of the parameters is reduced by assuming $`a_{11}=a_{22}`$, $`m_1=m_2`$ and $`\omega _1=\omega _2=\overline{\omega }`$, and $`N_1^{2\mathrm{d}}=N_2^{2\mathrm{d}}=N^{2\mathrm{d}}`$. For the typical experimental conditions $`(\omega _{},\omega _z)=2\pi (8,5)`$Hz and $`N10^6`$, the values of the intracomponent coupling constants are chosen as $`C_{11}=C_{22}=C=2000`$. In the numerical calculation, the equilibrium solutions are found by the norm-preserving imaginary time propagation of the time-dependent coupled GP equations, starting from arbitrary trial wave functions without vorticity. The propagation continues until the fluctuation in the norm of the wave function becomes smaller than 10-8.

With a fixed value $`C=2000`$, we investigate the equilibrium solutions for various values of two free parameters $`\delta C_{12}/C`$ and $`\overline{\mathrm{\Omega }}\mathrm{\Omega }/\overline{\omega }`$. The obtained vortex structure is summarized in the phase diagram of Fig. 1. Because the centrifugal potential overcomes the harmonic trapping potential for $`\overline{\mathrm{\Omega }}>1`$, the upper limit of the rotation frequency is set by $`\overline{\mathrm{\Omega }}=1`$. Previous studies have revealed that two non-rotating components with equal mass and equal intracomponent scattering length overlap completely for $`\delta <1`$, and phase-separate for $`\delta >1`$ Ho ; Ersy ; Tim . This criterion is reflected on the structure of vortex states, as explained below. In general, there are many metastable states in rotating two-component BECs; they are obtained by the numerical simulation starting from different wave functions. Here we discuss some characteristic features of vortex structures.

In the region $`\delta <1`$, two types of the regular vortex lattices are obtained as the equilibrium state. For $`\delta =0`$ where two components are not interacting, the formulation is equivalent to that of one component, thereby triangular vortex lattices are formed Tsubota . As $`\delta `$ increases, the positions of vortex cores in one component gradually shift from those of the other component and the triangular lattices are distorted. Eventually, the vortices in each component form a square lattice rather than a triangular one. The $`\delta `$ dependence on the stable region of a square lattice was studied in Ref. Mueller  in the high rotation limit $`\overline{\mathrm{\Omega }}1`$. As seen in Fig. 1, however, that stable region depends not only on $`\delta `$ but also sensitively on $`\overline{\mathrm{\Omega }}`$. An increase in rotation frequency also causes the transition from triangular to square lattices; in Figs. 2 (a) and (b), the two-dimensional density profiles of the condensates with $`\delta =0.7`$ are shown for $`\overline{\mathrm{\Omega }}=0.6`$ and $`\overline{\mathrm{\Omega }}=0.75`$. We find that the two vortex lattices are interlocked in such a manner that a peak in the density of one component is located at the density hole of the other, as shown in Fig. 2 (c). As a result, the total density $`\rho _T=|\mathrm{\Psi }_1|^2+|\mathrm{\Psi }_2|^2`$ obeys the Thomas-Fermi distribution applied to the overlapping two-component BECs with solid-body rotation $`\rho _T(r)=2\sqrt{\alpha /\pi }\alpha r^2`$ with $`\alpha =(1\overline{\mathrm{\Omega }}^2)/C(1+\delta )`$. We have confirmed that this feature holds for other parameter regimes.

Why is the square lattice stabilized in two-component BECs? According to the energy functional of this system, two components interact via the intercomponent interaction $`C_{12}|\mathrm{\Psi }_1|^2|\mathrm{\Psi }_2|^2`$; the velocity field in one component is independent of that of the other. Therefore, such a feature is determined only by the density distribution of the condensates which minimizes the interaction energy. We rewrite the interaction energy in terms of the total density $`\rho _T`$ and the “spin” variable $`S=|\mathrm{\Psi }_1|^2|\mathrm{\Psi }_2|^2`$ as $`E_{\mathrm{int}}=𝑑𝐫C[(1+\delta )\rho _T^2+(1\delta )S^2]/4`$. A larger $`\delta `$ makes a smoother total density $`\rho _T`$ more favorable in order to reduce the first term of $`E_{\mathrm{int}}`$. This results in the shift of the positions of vortex cores of each component. Then, the last term of $`E_{\mathrm{int}}`$ can be interpreted effectively as the interaction energy between spins in the Ising model atten2 . Spin-up components correspond to the density peaks of $`\mathrm{\Psi }_1`$ at the vortex cores of $`\mathrm{\Psi }_2`$, and vice versa for spin-down components. When the coefficient $`1\delta `$ is positive, the interaction between spins is anti-ferromagnetic, which makes a square lattice stabilized (because the triangular lattice could be frustrated). Anti-ferromagnetic nature is made more pronounced as the interlocking of vortex lattices becomes stronger with $`\delta `$ and the number of vortices increase with $`\overline{\mathrm{\Omega }}`$.

As $`\delta `$ exceeds unity, the system enters a ferromagnetic phase. Then, the condensates undergo phase separation to spontaneously form domains having the same spin component atten2 . Concurrently, vortices of the same component begin to overlap at $`\delta 1`$. Figure 3 shows typical solutions at $`\delta =1`$. In Fig. 3(a), vortices of each component are overlapped in lines, and each condensate density forms a stripe pattern Engels . While this structure can also be derived if one assumes a perfect lattice Mueller , we find that its energy is nearly degenerate with that of Fig. 3(b). For the same parameters, we can obtain another equilibrium state called “double-core vortex lattice” in Fig. 1, where a vortex lattice of component 2 is made by pairs of vortices with the same circulation; vortices in component 1 surround those pairs. Therefore, various metastable structures will appear in this parameter region.

For strongly phase-separated region $`\delta >1`$, the density peaks of the same spin component, at which the other-component vortices are located, merge further, resulting in the formation of vortex sheets; a typical example is shown in Fig. 4(a). Singly quantized vortices line up in sheets, and the sheets of component 1 and 2 are interwoven alternately from the center to the outward. Figure 4(b) shows that each superfluid velocity $`v_i`$ $`(i=1,2)`$ jumps at the vortex sheet, following approximately the velocity of solid-body rotation. In the region $`\delta >1`$, even though the value of $`\delta `$ is changed, the total density is fixed by the Thomas-Fermi distribution with $`\delta =1`$.

The stationary vortex sheet has been observed in rotating superfluid 3He-A sheethe , where vortices are bounded to a topologically stable domain wall soliton across which the unit vector representing the direction of the orbital angular momentum of the Cooper-pair faces oppositely. As discussed in superfluid 3He-A, the equilibrium distance between the sheets of component 1 and component 2 is determined by the competition between the surface tension $`\sigma `$ of the domain wall and the kinetic energy of the superflow sheethe . In order to estimate that distance in two-component BECs, we consider a simple model shown in Fig. 4(c). In this model, each velocity is assumed to be constant between the vortex sheets of the same component; the value of $`v_i`$ increases by $`2\overline{\mathrm{\Omega }}b`$ across every sheet. The total density $`\rho _T=|\mathrm{\Psi }_1|^2+|\mathrm{\Psi }_2|^2`$ is constant, and the domain boundary with the penetration depth $`\mathrm{\Lambda }(<b)`$ is approximated by the linear profile. We calculate the free energy $`F=E\mathrm{\Omega }L_z`$ in the range $`0<r<2b`$, the sheet distance $`b`$ being determined so as to minimize $`F`$ per unit area. Firstly, the penetration depth $`\mathrm{\Lambda }`$ is determined by minimizing the surface tension $`\sigma `$ of a single domain wall with respect to $`\mathrm{\Lambda }`$ Tim , which is the sum of the quantum pressure energy $`E_{\mathrm{qp}}/2\pi b(\rho _T/2\mathrm{\Lambda })`$ and the interaction energy in the overlap region $`E_{\mathrm{int}}/2\pi b=(\delta 1)\mathrm{\Lambda }\rho _T^2/12`$; we thus obtain $`\mathrm{\Lambda }_p=\sqrt{6/C(\delta 1)\rho _T}`$. Then, the surface tension is written as $`\sigma =\rho _T^{3/2}\sqrt{C(\delta 1)/6}`$. Secondly, the flow energy per area (in rotating frame) is given by $`(1/4\pi b^2)_0^{2b}d^2r_i\rho _i(v_i\overline{\mathrm{\Omega }}r)^2/2=29\rho _T\overline{\mathrm{\Omega }}^2b^2/768`$, where $`b\mathrm{\Lambda }_p`$ is assumed. Thus, the free energy per unit area is written as

$$\frac{F}{4\pi b^2}=\frac{29\rho _T\overline{\mathrm{\Omega }}^2b^2}{768}+\frac{2\sigma }{b}+\frac{C\rho _T^2}{2}.$$
(1)
Minimizing this energy with respect to $`b`$, one obtains $`b=(768\sigma /29\overline{\mathrm{\Omega }}^2\rho _T)^{1/3}=(768/29\overline{\mathrm{\Omega }}^2\mathrm{\Lambda }_p)^{1/3}`$. By using the Thomas-Fermi density at $`r=0`$ with $`\delta =1`$ as the value of $`\rho _T`$, the sheet spacing $`b(\delta 1)^{1/6}(1\overline{\mathrm{\Omega }}^2)^{1/12}/\overline{\mathrm{\Omega }}^{2/3}`$ is consistent with that of numerical solutions; for example, for parameters used in Fig. 4 we obtain $`b=3.09b_{\mathrm{ho}}`$. The vortex sheet is expected in the region $`b>\mathrm{\Lambda }_p`$ and $`2b<R_{\mathrm{TF}}=\sqrt{2/\sqrt{\pi \alpha }}`$, as shown in the shaded region in Fig. 1. When $`2b>R_{\mathrm{TF}}`$, i.e., $`b`$ becomes comparable with the condensate size, the clear structure of sheets vanishes.

Finally, we comment on what structure could be observed in actual experimental conditions. The two-component BECs realized in JILA is a mixture of the states $`|1,1`$ and $`|2,1`$ of 87Rb Hall . This mixture has the scattering lengths which have the same order, $`a_{11}:a_{22}:a_{12}=1:0.94:0.97`$, i.e., $`\delta 1`$ for equal particle numbers. Therefore, lattices with partially overlapping vortices of the same component, as shown in Fig. 3, are expected to be observed. The MIT group has observed the phase separation of a mixture of 23Na BECs with $`|1,0`$ and $`|1,1`$ state Miesner , which has $`\delta 1`$ ($`a_{11}:a_{22}:a_{12}=1:1.035:1.035`$). In addition, a mixture of 41K and 87Rb BECs reported recently Modugno  lies deeply in a phase-separate region. The vortex sheets should therefore be observed at high rotation frequencies.

In conclusion, we reveal a rich phase diagram of vortex states in rotating two-component BECs, even in the restricted parameter space $`(\delta ,\overline{\mathrm{\Omega }})`$. The structure of a vortex lattice is shown to depend sensitively on the rotation frequency. In particular, we find new phenomena such as double-core vortex lattices and interwoven serpentine vortex sheets. Use of different atomic masses and intracomponent interactions will realize a coexistence system of vortices with different vortex-core sizes; such a situation may change the lattice structure shown in this Letter. We plan to study a more detailed phase diagram and effects of the internal Josephson coupling on the vortex states.

K.K. and M.T. acknowledge support by a Grant-in-Aid for Scientific Reserach (Grants No. 1505955 and No. 15341022) by the Japan Society for the Promotion of Science. M.U. acknowledges support by a Grant-in-Aid for Scientific Research (Grant No.14740181) by the Ministry of Education, Culture, Sports, Science and Technology of Japan, and a CREST program by Japan Science and Technology Corporation (JST).


format: text

text:
# An Application of Feynman-Kleinert Approximants to the Massive Schwinger Model on a Lattice

## I INTRODUCTION

The analytic continuation of a function, whose perturbation series expansion about a given point is known, is an important and long-standing problem. Some standard methods are the use of Padé approximants, or integrated differential approximants, as reviewed for instance by Guttmann. In some circumstances, these techniques can work extremely well; but in others they fail. A well known difficult case is that of a function whose series expansion about the origin in a variable $`x`$, say, is known up to a given order, and we want to extrapolate to $`x\mathrm{}`$. If there is a cut or an essential singularity at infinity, the standard methods do not work well. This is a common situation in lattice gauge theory, for instance.

In recent years, a new method has been developed, involving “Feynman-Kleinert approximants”. This is a variational technique, arising originally from a Gaussian approximation to the path integral. It has been used with great success on a number of problems, such as extrapolation of the weak-coupling series for the anharmonic oscillator, where the asymptotic strong-coupling behaviour was obtained with extraordinary accuracy, up to 20 significant digits. For a discussion and further references, see the recent paper by Kleinert.

Here we attempt to apply the new method to a lattice gauge theory, namely the Schwinger model in (1+1) dimensions. In the series approach to lattice gauge theory, one generates strong-coupling series expansions for the physical quantities of interest, and then tries to extrapolate to the continuum, weak-coupling limit. It would be an important breakthrough if a better technique could be found for performing this extrapolation. In the report we set out to test whether Feynman-Kleinert approximants can fulfill this need. Unfortunately, we find that the technique is not very successful in this context.

Another common problem is that of interpolation of a function, where a number of series coefficients are known of both the weak-coupling and the strong-coupling expansions. An example of this arises for the low-lying energy eigenvalues in the continuum Schwinger model, as functions of $`m/g`$, where $`m`$ is the fermion mass and $`g`$ is the electric coupling. We show that Feynman-Kleinert approximants can provide an effective interpolation between these two limits.

## II THEORY

For reference, let us briefly recapitulate the algorithm set out by Kleinert using his notation. Consider a quantity $`E`$ whose perturbation series expansion in some parameter $`\alpha `$ is known to order $`N`$:

$$E_N=\underset{n=0}{\overset{N}{}}a_n\alpha ^n.$$
(1)
Write

$$E_N^\omega =\omega ^p\underset{n=0}{\overset{N}{}}a_n\left(\frac{\alpha }{\omega ^q}\right)^n,$$
(2)
where $`\omega `$ is an auxiliary parameter whose value will eventually be set equal to 1, and $`p,q`$ are indices which will be chosen later to suit the problem at hand. Now set $`\omega \sqrt{\mathrm{\Omega }^2+\omega ^2\mathrm{\Omega }^2}`$, and re-expand $`E_N^\omega `$ in powers of $`\lambda =1\omega ^2/\mathrm{\Omega }^2`$, treating $`\lambda `$ as a quantity of order $`\alpha `$, and truncating the re-expanded series at order $`N`$. The result is a function

$$W_N(\alpha ,\mathrm{\Omega })=\mathrm{\Omega }^p\underset{n=0}{\overset{N}{}}a_nf_n(\mathrm{\Omega })\left(\frac{\alpha }{\mathrm{\Omega }^q}\right)^n,$$
(3)
where

$$f_n(\mathrm{\Omega })=\underset{j=0}{\overset{Nn}{}}\left(\begin{array}{c}(pqn)/2\\ j\end{array}\right)(1)^j\left(1\frac{\omega ^2}{\mathrm{\Omega }^2}\right)^j.$$
(4)
Forming the first and second derivatives of $`W_N(\alpha ,\mathrm{\Omega })`$ with respect to $`\mathrm{\Omega }`$ (and setting $`\omega =1`$), we find the positions of the turning points. The smallest among these is denoted $`\mathrm{\Omega }_N`$. The resulting $`W_N(\alpha )W_N(\alpha ,\mathrm{\Omega }_N)`$ gives the desired approximation to the function $`E`$ at finite $`\alpha `$, an extrapolation of the function $`E_N`$ to finite $`\alpha `$.

Now consider the limit $`\alpha \mathrm{}`$. At large $`\alpha `$, $`\mathrm{\Omega }_N`$ will scale with $`\alpha `$ like

$$\mathrm{\Omega }_N\alpha ^{1/q}c,$$
(5)
so that

$$W_N(\alpha ,\mathrm{\Omega }_N)\alpha ^{p/q}.$$
(6)
The full ‘strong-coupling’ (large $`\alpha `$) expression can now be obtained by writing

$$W_N(\alpha ,\mathrm{\Omega })=\mathrm{\Omega }^pw_N(\widehat{\alpha },\omega ^2/\mathrm{\Omega }^2),$$
(7)
with $`\widehat{\alpha }=\alpha /\mathrm{\Omega }^q`$, and expanding $`w_N`$ in powers of $`\omega ^2/\mathrm{\Omega }^2`$, which for $`\alpha \mathrm{}`$ behaves like $`(1/c^2)(\alpha /\omega ^q)^{2/q}`$. The result is

$$W_n(\alpha )=\alpha ^{p/q}\left[b_0(c)+b_1(c)\left(\frac{\alpha }{\omega ^q}\right)^{2/q}+b_2(c)\left(\frac{\alpha }{\omega ^q}\right)^{4/q}+\mathrm{}\right]$$
(8)
with

$`b_n(c)`$ $`=`$ $`{\displaystyle \frac{1}{n!}}w_N^{(n)}(\widehat{\alpha },0)\widehat{\alpha }^{(2np)/q}|_{\widehat{\alpha }=1/c^q}`$ (9)
$`=`$ $`{\displaystyle \underset{l=0}{\overset{N}{}}}a_l{\displaystyle \underset{j=n}{\overset{Nl}{}}}\left(\begin{array}{c}(plq)/2\\ j\end{array}\right)\left(\begin{array}{c}j\\ n\end{array}\right)(1)^{jn}c^{plq2n}.`$ (14)
The indices $`p`$ and $`q`$ in the expansion (2) are now chosen so as to give the correct asymptotic behaviour (i.e. the correct powers of $`\alpha `$) in the strong-coupling expansion (8), which we assume is known a priori. The leading coefficient $`c`$ in (5) is found by searching for the extrema of the leading coefficient $`b_0(c)`$ as a function of $`c`$ and choosing the smallest of them.

Finally, we have to account for the fact that $`\mathrm{\Omega }_N`$ will have corrections to the leading behaviour $`\alpha ^{1/q}c`$. We can allow for this by replacing $`c`$ by

$$c(\alpha )=c_0+c_1\left(\frac{\alpha }{\omega ^q}\right)^{2/q}+c_2\left(\frac{\alpha }{\omega ^q}\right)^{4/q}+\mathrm{}$$
(15)
requiring a re-expansion of the coefficients $`b_n(c)`$ in (14). The expansion coefficients $`c_n`$ are determined by looking for the turning points of $`b_{2n}(c)`$, successively.

The equations

$$\frac{db_{2n}(c)}{dc_n}=0,$$
(16)
together with knowledge of the coefficients $`\{a_n,n=0,\mathrm{},N\}`$, can in principle be used to determine a set of parameters $`\{c_n,n=0,\mathrm{},N\}`$, and hence estimate the leading strong-coupling coefficients $`b_n(c)`$. This provides estimates of the asymptotic behaviour of the function $`E`$. In the case of the strong-coupling anharmonic oscillator, it has been shown that the resulting estimates converge exponentially fast, and have been used to estimate the energy eigenvalue in the strong-coupling limit accurate to 20 decimal places!

Alternatively, any knowledge of the strong-coupling coefficients $`b_0,b_1,\mathrm{}`$ can be fed back into the set of equations (14) and (16) to obtain estimates of further weak-coupling coefficients $`a_{N+1},a_{N+2},\mathrm{}`$, and thus carry the extrapolation $`W_N(\alpha )`$ to higher orders. This is the basis of the interpolation algorithm between weak and strong-coupling for the function $`E`$.

## III RESULTS

### A Extrapolation to the Continuum Limit

We have applied the Feynman-Kleinert approximation method to the lattice Schwinger model, as formulated by Banks et al. . The rescaled lattice Hamiltonian is

$$W=\frac{2}{ag^2}H=W_0+xV,$$
(17)
where

$`W_0={\displaystyle \underset{n}{}}L^2(n)+\mu {\displaystyle \underset{n}{}}(1)^n\varphi ^{}(n)\varphi (n),`$ (18)
$`V=i{\displaystyle \underset{n}{}}\left[\varphi ^{}(n)e^{i\theta (n)}\varphi (n+1)\text{h.c.}\right],`$ (19)
and

$$\mu =\frac{2m}{g^2a},x=\frac{1}{g^2a^2}.$$
(20)
Here $`m`$ is the electron mass, $`g`$ is the electric coupling constant, and $`a`$ is the lattice spacing. The field $`\varphi (n)`$ is the lattice fermion field defined on sites $`n`$,

$$\{\varphi (n),\varphi ^{}(m)\}=\delta _{nm},$$
(21)
while $`L(n)`$ is the electric field on link $`n`$, which takes only integer eigenvalues.

Expansions of the energy eigenvalues of the operator $`W`$ have been calculated to high orders in the perturbation parameter $`x`$, using efficient linked cluster methods. Our object now is to see whether these series can be accurately extrapolated to the continuum limit, which corresponds to $`a0`$ or $`x\mathrm{}`$, using the Feynman-Kleinert approach.

#### 1 Ground State Energy

For the ground state energy per site $`\omega _0/M`$, where $`M`$ is the number of states, the series expansion is of the form

$$\omega _0/M=\underset{n}{}a_n\left(x^2\right)^n,$$
(22)
where $`a_0=0`$; and the asymptotic behaviour as $`x\mathrm{}`$ is known

$$\omega _0/M\frac{2x}{\pi }.$$
(23)
The coefficients $`a_n`$ have been calculated up to order $`n=15`$ . In this case, we can apply the formalism of section II, taking $`E=\omega _0/M`$, $`\alpha =x^2`$, $`p=2`$, $`q=4`$, which implies the asymptotic behaviour as $`x\mathrm{}`$

$$\omega _0/Mb_0x+b_1+O(x^1).$$
(24)
Fig. 1 shows the resulting approximants to the ground-state energy at $`m/g=0`$ compared with numerical results obtained by other methods, and the known continuum limit, Eq. (23). It can be seen that the approximants do seem to converge to the correct limit, but only rather slowly.

#### 2 Energy Gap

The energy gap to the lowest-lying positronium state has an expansion

$$\omega _1\omega _0=\underset{n}{}a_n\left(x^2\right)^n$$
(25)
and the expected asymptotic behaviour as $`x\mathrm{}`$ is

$$\omega _1\omega _0b_0x^{1/2}.$$
(26)
The coefficients $`a_n`$ for this series have been calculated up to order $`n=14`$ . Here, we take $`\alpha =x^2`$, $`p=2`$, $`q=8`$, assuming the asymptotic behaviour as $`x\mathrm{}`$

$$\omega _1\omega _0b_0x^{1/2}+b_1+O(x^{1/2})$$
(27)
which accords with previous studies. In this case, however, the Feynman-Kleinert approximants do not converge at all at large $`x`$. Fig. 2 shows the dependence of the estimated value of $`b_0`$ as a function of the parameter $`c_0`$, for the case $`m/g=0`$, at various orders $`N`$. It can be seen that as $`c_0`$ increases, the estimate of $`b_0`$ comes in from infinity, oscillates around the correct value (marked by a solid line ), and then slowly drifts away again. Unfortunately, the oscillations increase with order $`N`$, and the values at the lowest turning point correspondingly diverge away from the correct value as $`N`$ increases. The same is true of the values at the lowest point of inflection. It seems that the Feynman-Kleinert approximant method fails in this case.

### B Interpolation of Continuum Series

Another problem connected with series expansions for the Schwinger model arises in connection with the continuum field theory, rather than the lattice model. A few coefficients are known of the series expansions for the lowest-lying excited states of the model, the “vector” and “scalar” positronium states, in both the strong-coupling and weak-coupling regimes. At strong coupling, $`m/g0`$, the first three coefficients have been computed by Adam, using the bosonic form of the theory:

$$\frac{E}{g}=\left(\frac{M2m}{g}\right)=\underset{n}{}a_n\left(m/g\right)^n.$$
(28)
At weak coupling, $`m/g\mathrm{}`$, the first three coefficients have been calculated by Sriganesh et al. , following a treatment of Coleman :

$$\frac{E}{g}\left(\frac{g}{m}\right)^{1/3}\underset{n}{}b_n\left[\left(\frac{g}{m}\right)^{2/3}\right]^n.$$
(29)
The problem now is, how can one compute a sensible interpolation between these two series to obtain accurate estimates of the energy eigenvalues at finite $`m/g`$? One possibility is to use “2-point” Padé approximants: but they involve only integer powers of the variable, $`(m/g)`$ in this case. By fitting an appropriate power of the function $`E/g`$ one can mimic the leading behaviour as $`m/g\mathrm{}`$, but cannot reproduce the sub-leading powers. Furthermore, experience shows that 2-point Padé approximants are prone to exhibit spurious ‘wiggles’ at intermediate couplings between the two limits. Here we set out to apply Feynman-Kleinert approximants to the problem.

Setting $`\alpha =m/g`$, $`p=1`$, $`q=3`$, we have computed Feynman-Kleinert approximants for the energy eigenvalues using known values for the coefficients $`\{a_n,b_n,n=0,N1\}`$ for $`N=1,2,3`$. Now knowing $`2N`$ coefficients in total allows us to estimate the next $`N`$ coefficients for either the $`a_n`$ or $`b_n`$ series. Table I shows the ‘predicted’ values versus the known values for $`N=1`$. For the vector state with $`N=1`$, the predicted values are in quite good agreement with the known values, to within a few percent. At $`N=2`$ (see Table II), however, a large discrepancy becomes evident between the predicted and ‘known’ values for the coefficient $`b_2`$ \- even the sign is different. This prompted us to re-examine the calculation of the weak-coupling coefficients $`b_n`$ by Sriganesh et al.; and indeed we discovered some errors in those calculations. A corrected version of the calculations is given in the Appendix; the corrected series coefficients are also shown in Table II.

For the scalar state, the ‘predicted’ series coefficients disagree with the known ones even at $`N=1`$ \- but this is easily understood, in that the simplest interpolation between the zeroth order series at either end would have no ‘hump’ in the middle such as is shown in Fig. 4. The predicted values at $`N=2`$ agree qualitatively with the (corrected) known values.

The Feynman-Kleinert approximants were found to give smooth and apparently accurate interpolations of the energy eigenvalues at finite $`m/g`$, up to order $`N=2`$. At order $`N=3`$, no solution was found to the system of simultaneous equations - we have not explored in detail why this occurred. Figures 3 and 4 show the results for the vector and scalar state interpolations. For the vector state energies, our interpolation is compared to data obtained from the density matrix renormalization group (DMRG) method , the “fast moving frame” estimates of Kröger and Scheu , and the “renormal-ordered” mass perturbation series results of Adam . The scalar state interpolations are compared to the finite lattice calculations of Sriganesh et al. , and again results of Kröger and Scheu, and Adam. Both plots show excellent agreement between the interpolation and the numerical data, within errors. We also plot the second order strong and weak coupling series used to generate the interpolation on the same plots, which diverge away for large and small masses respectively. The “renormal-ordered” expansion  about $`m/g=0`$, extends the region of validity of the perturbation theory, and agrees fairly well with our interpolation for small masses.

The interpolation performs particularly well for the scalar state, where for small $`m/g`$, the original strong coupling series disagrees significantly with the numerical data. The interpolation on the other hand agrees extremely well with existing results, despite the large gap that must be bridged between the two series. Tables III and IV list numerical values for specific values of $`m/g`$ for the interpolation compared to previous works. We estimate the error on our interpolations based on the difference between our quoted $`N=2`$ calculation and the corresponding $`N=1`$ calculation. Error estimates have probably been overestimated somewhat using this method, as can be seen from the good agreement with the numerical data in Figure 4. We have not however found a better way of estimating these errors.

## IV CONCLUSIONS

The Feynman-Kleinert approximants have earned a mixed scorecard in the applications we have discussed. For the extrapolation of lattice strong-coupling series to the weak-coupling continuum limit, the approximants have basically failed. They converged for the ground-state eigenvalue, but only slowly; and for the first excited state eigenvalue, they did not converge at all. This is probably due to the presence of exponentially decaying terms corresponding to essential singularities in the weak-coupling limit, which are not accounted for in the Feynman-Kleinert approach. The upshot is, unfortunately, that this approach does not seem to have any advantage over previous techniques, which used Padé or integral approximants to extrapolate to intermediate couplings, and then ‘matched’ the results onto a weak coupling form to reach the continuum limit.

For the continuum field theory, however, the Feynman-Kleinert approximants have proved useful in providing a sensible and apparently accurate interpolation between the known weak coupling and strong coupling series. The fact that we were able to pick up a mistake in the series coefficients by this approach adds confidence in the technique. No doubt many applications of a similar nature are possible.

###### Acknowledgements.

We are grateful to Dr. C. Adam for very useful correspondence on this topic. This work forms part of a research project supported by the Australian Research Council.

## A WEAK COUPLING EXPANSION OF THE MASSIVE SCHWINGER MODEL

Our starting point is the Hamiltonian as derived by Coleman in the two-particle subspace (see Eq. (4.18) of Ref. ). We set the background field $`\theta =0`$. Expanding this Hamiltonian for large mass $`m`$, and keeping terms to order $`O(m^3)`$ we obtain a Schrödinger equation

$$(H_0+H_1)\psi (x)=E\psi (x),$$
(A1)
where

$`H_0`$ $`=`$ $`{\displaystyle \frac{1}{m}}{\displaystyle \frac{^2}{x^2}}+{\displaystyle \frac{g^2}{2}}|x|`$ (A2)
$`H_1`$ $`=`$ $`{\displaystyle \frac{g^2}{m\pi }}{\displaystyle \frac{1}{4m^3}}{\displaystyle \frac{^4}{x^4}}+{\displaystyle \frac{g^2}{4m^2}}\delta (x){\displaystyle \frac{g^2}{2\pi m^3}}{\displaystyle \frac{^2}{x^2}}.`$ (A3)
Rescale to dimensionless variables as follows:

$$x=\left(\frac{2}{mg^2}\right)^{1/3}zE=\left(\frac{g^4}{4m}\right)^{1/3}ϵ.$$
(A4)
Then (A1) becomes

$$(H_0^{}+H_1^{})\psi (z)=ϵ\psi (z),$$
(A5)
where

$`H_0^{}`$ $`=`$ $`{\displaystyle \frac{^2}{z^2}}+|z|`$ (A6)
$`H_1^{}`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}\left({\displaystyle \frac{2g}{m}}\right)^{2/3}\left({\displaystyle \frac{g}{4m}}\right)^{4/3}[{\displaystyle \frac{^4}{z^4}}2\delta (z)].`$ (A7)
The last term in (A2) has now been dropped as it is now clear it is of order $`(g/m)^2`$. The solution of the leading-order Schrödinger equation

$$\left(\frac{d^2}{dz^2}|z|\right)\psi (z)=ϵ\psi (z),$$
(A8)
was discussed by Hamer. The solution of (A8) is a symmetric or antisymmetric Airy function which obeys the condition

$`\text{Ai}^{}(ϵ_n)=0`$ (symmetric) (A9)
$`\text{Ai}(ϵ_n)=0`$ $`\text{(antisymmetric)}.`$ (A10)
The lowest symmetric (antisymmetric) wavefunction gives the energy for the vector (scalar) state.

The higher order terms contained in $`H_1^{}`$ may be taken into account by calculating the correction term $`\mathrm{\Delta }ϵ=\psi |H_1^{}|\psi `$ using numerical integration (the results are shown in Table V). Rescaling back to our original variables gives us

$$\frac{E_1}{g}=0.6418\left(\frac{g}{m}\right)^{1/3}\frac{1}{\pi }\left(\frac{g}{m}\right)+0.1547\left(\frac{g}{m}\right)^{5/3},$$
(A11)
for the vector state binding energy, and

$$\frac{E_2}{g}=1.473\left(\frac{g}{m}\right)^{1/3}\frac{1}{\pi }\left(\frac{g}{m}\right)0.1093\left(\frac{g}{m}\right)^{5/3}.$$
(A12)
for the scalar state binding energy.


text:
# Precision Cosmology? Not Just Yet… Sarah L. Bridle, Ofer Lahav, Jeremiah P. Ostriker, Paul J. Steinhardt 11footnote 1S. L. Bridle J. P. Ostriker and O. Lahav are at the Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK. P. J. Steinhardt is in the Department of Physics, Princeton University, Princeton, New Jersey 08544, USA.

The recent announcement by the WMAP satellite team of their landmark measurements of the cosmic microwave background (CMB) anisotropy *(1-3)* has convincingly confirmed important aspects of the current standard cosmological model. The results show with high precision that space is flat (rather than curved) and that most of the energy in the universe today is “dark energy”, which is gravitationally self-repulsive and accelerates the expansion of the universe. The evidence is independent of supernovae results *(4,5)*.

The measurements strongly indicate that the amplitudes of spatial variations in density and temperature that seeded the formation of galaxies were roughly independent of length scale, adiabatic (all forms of energy have the same spatial variation), and followed a Gaussian distribution — just as predicted by the standard Big Bang inflationary model. WMAP heralds a new age of precision cosmology with careful error analysis, tightly constraining many key parameters *(6)*. For example, the lifetime of the universe has been determined to be $`13,400\pm 300`$ million years *(6)*. Furthermore, WMAP’s new measurement of the CMB polarization as a function of angular scale shows that the epoch of cosmic reionization — associated with the formation of the first stars — had already occurred when the universe was several hundred million years old.

At the same time we celebrate this triumph, it is important to recognize that important issues remain. For example, it is not yet clear whether the spectrum of temperature fluctuations is truly consistent with inflation. The spectrum is roughly scale-invariant, but there are hints of peculiarities, and a key inflationary prediction — the presence of gravitational wave effects — has not yet been observed.

We also do not know whether dark energy is due to an unchanging, uniform, and inert “vacuum energy” (also known as a cosmological constant) or a dynamic cosmic field that changes with time and varies across space (known as quintessence). “Dark matter”, which is gravitationally self-attractive, also remains mysterious: We do not yet know its nature, nor are we certain about its density or the amplitude of the initial ripples in its distribution.

Today’s standard theoretical paradigm is the inflationary Big Bang model. According to this picture, the universe began in a a state of nearly infinite temperature and density and almost immediately entered a phase of rapid, accelerated expansion (“inflation”). This expansion smoothed out the distribution of energy, flattened initial warp or curvature in space, and created tiny variations in density. To transform these density variations into the gravitationally collapsed, complex structures we see today, it is essential that there be “dark matter” as well as ordinary (baryonic) matter. Finally, we need dark energy to account for the measured total energy density and to explain the current cosmic acceleration.

Some of the WMAP results – the flatness of space, the near scale-invariance, adiabaticity, and Gaussian distribution of the density perturbations *(7)*, the density of baryons, the age of the universe, and perhaps the early formation of the first stars — are based on WMAP alone and are consistent with the standard model. Because the CMB is a direct image of the early universe and its interpretation entails simple, well-understood physical principles, these results are robust.

On the other hand, some important issues can only be addressed by combining WMAP data with other cosmological measurements. These conclusions should be viewed more cautiously because the result depends sensitively on the choice of additional data.

For example, by combining data, a significant deviation from a perfectly scale-invariant $`(n=1)`$ spectrum was found *(8)*. According to the best-fit WMAP combined analysis *(8)*, $`n`$ runs from 1.1 on the largest scales to $`<0.9`$ on the smallest scales probed, a deviation that disagrees with the simplest and most natural inflationary models *(9)*. These results cast a pall over the inflationary paradigm at the same time that many of its other important predictions are confirmed. However, it is important to note that WMAP data alone are not inconsistent with the simplest inflationary models, the fit being $`n=0.99\pm 0.04`$. The inconsistency only arises when the data are combined with two degree field galaxy redshift survey and quasar absorption line measurements.

Setting aside this issue, we next consider whether the WMAP data are uniquely consistent with the standard inflationary Big Bang picture. The answer is no, as the WMAP team has itself indicated. There remains room for radical alternatives. An example is the cyclic picture, in which the universe undergoes a periodic sequence of cycles *(10)*: Expansion from a hot big bang is followed by contraction in a “big crunch” and reemergence in a big bang, and the key events that shaped the large-scale structure of the universe occur before the “bang”, a cycle ago. This model offers a very different view of cosmic history, yet it fits all current observations (including the new WMAP results) at least as well as the inflationary picture.

Even if the standard picture is proven to be correct, the model is incomplete. The initial conditions that led to inflation and the identity of the “inflaton” field (the cosmic field that causes inflationary expansion) remain unknown and the nature of dark matter is unsettled.

An important uncertainty is the nature of the gravitationally self-repulsive dark energy. Whether it consists of vacuum energy or quintessence depends on $`w`$, the ratio of pressure to energy density. The WMAP combined analysis concluded that the best fit corresponds to vacuum energy. But their own closer examination of WMAP results at large angular scales shows that the enhanced fluctuations expected for vacuum energy are missing. This discrepancy could be a first sign that dark energy is quintessential *(11)*. The combined analysis washes out the effect, but this could be hiding a very important hint about the true nature of dark energy and confounding measurements of parameters.

WMAP has produced impressive constraints on many fundamental cosmic parameters \[see *(6)* for extensive tables\], but there remain open frontiers. Perhaps most important is the uncertainty in the density of dark matter, $`\mathrm{\Omega }_\mathrm{m}`$, and the amplitude of density fluctuations labeled by $`\sigma _8`$, as represented in the figures. These parameters are critical because they describe the amount and distribution of the matter that clusters to form all of the structure in the universe.

The WMAP data alone only determine some combination of $`\mathrm{\Omega }_\mathrm{m}`$ and $`\sigma _8`$, represented by the narrow pink contour in the first figure, if one takes into account the uncertainties in just three other cosmological parameters but assumes the rest are perfectly known (“strong prior”). Even with these over optimistic assumptions, the range of $`\mathrm{\Omega }_\mathrm{m}`$ and $`\sigma _8`$ is large. As we progressively relax the assumptions by including empirical uncertainties, the range of degeneracy balloons until we get to the more realistic purple contour (the “weak prior”).

To reduce the uncertainty, the WMAP team has presented results obtained by combining with the two degree field galaxy redshift survey and quasar absorption line measurements. We agree with the mathematical conclusions obtained by the WMAP team when combining these data sets. On the other hand, it is worth considering the broader range of available data shown in the second figure, which reveals three different directions of degeneracy: vertical strips from methods that measure the matter density alone; constraints on the combination $`\sigma _8\mathrm{\Omega }_\mathrm{m}^{0.6}`$, obtained by measuring the number of galaxy clusters observed today and the apparent distortion of galaxy images due to the bending of light by dark matter; and a roughly orthogonal constraint from the WMAP data. Here we see that the high likelihood regions (solid lines) do not all overlap well, suggesting a problem with one or more of the measurements, or their interpretation, or, more interestingly, that the underlying model may be wrong.

Perhaps adding only select measurements to WMAP data will prove to be the correct strategy. On the other hand, given the issues raised in the second figure, it may be that the real uncertainty is much greater. We will have to wait for forthcoming observations to substantially reduce the current uncertainty.

Thus, even after the historic WMAP breakthrough, there remain unresolved issues, and so there is plenty of room for surprises. Only 5 years ago, breakthroughs in technology and astronomical technique led to the discovery that the expansion of the universe is accelerating. The future holds promise of even greater technological advances that will uncover further cosmological surprises.

1. C. Bennett *et al.* astro-ph/0302207
2. G. Hinshaw *et al.* astro-ph/0302217
3. A. Kogut *et al.* astro-ph/030221
4. S. Perlmutter *et al.*, ApJ 517, 565 (1999)
5. A.G. Riess, *et al.*, AJ, 116, 1009 (1998)
6. D. Spergel *et al.* astro-ph/0302209
7. E. Komatsu *et al.* astro-ph/0302223
8. H. Peiris *et al.* astro-ph/0302225
9. J. Khoury, P. Steinhardt, N. Turok astro-ph/0302012
10. P. Steinhardt, N. Turok Science 296, 1436 (2002)
11. S. DeDeo, R. Caldwell, P. Steinhardt astro-ph/0301284
12. L. Verde *et al.* astro-ph/0302218
13. C. Kuo *et al.* astro-ph/0212289
14. T. Pearson *et al.* astro-ph/0205388
15. K. Grainge *et al.* astro-ph/0212495
16. A. Lewis, S. Bridle PRD 66, 103511 (2002)
17. N. Bahcall, R. Cen, R.Davé, J.Ostriker, Q.Yu ApJ 541, 1 (2000)
18. S. Allen, R. Schmidt, A. Fabian, H. Ebeling astro-ph/0208394
19. E. Pierpaoli, S. Borgani, D. Scott, M. White astro-ph/0210567
20. H. Hoekstra, H. Yee, M. Gladders ApJ 577, 595 (2002)
21. M. Jarvis *et al.* astro-ph/0210604
22. We thank Antony Lewis for making the CosmoMC code publicly available and for helping with the WMAP constraints. Computations were performed at the UK National Cosmology Supercomputer Center funded by PPARC, HEFCE and Silicon Graphics / Cray Research. Cluster counts constraints were kindly provided by S. Allen. We thank the Leverhulme Quantitative Cosmology group in Cambridge for helpful discussions.


text:
# Averages and critical exponents in type-III intermittent chaos

## I Introduction

Intermittent chaos is the phenomenon shown by systems exhibiting long sequences of periodic-like behaviour, the laminar phases, separated by comparatively short chaotic eruptions. Intermittent chaotic systems have been extensively studied since the original proposals by Pomeau and Manneville pomeau80 ; berge84 ; manneville90 ; kim98  classifying type-I, II, and III instabilities when the Floquet multipliers of a map crosses the unity circle. Although other mechanisms may occur leading to intermittency these three cases are the most simple and most frequently encountered in low dimensional systems. The class of intermittent chaos studied by Pomeau and Manneville includes the tangent bifurcation, leading to intermittency of type I, when the Floquet’s Multiplier for the associated map crosses the circle of complex numbers with unitary norm through $`+1`$; the Hopf bifurcation, leading to type-II intermittency, which appears as two complex eigenvalues of the Floquet’s Matrix cross the unitary circle off the real axis; and the sub-critical period doubling, leading to type-III intermittency, whose critical Floquet’s Multiplier is $`1`$.

Many experimental evidences for these intermittent behavior have appeared in the literature. The type-III has been reported for electronic nonlinear devices jeffries82 ; fukushima88 ; ono95 , lasers tang91 , and biological tissues griffith97 . A signature for this intermittency is given by the critical exponent describing the dependence of average length of nearly periodic phases, that is, laminar phases, with the control parameter. Histograms of number of laminar phases longer than a given duration are related to the normal form nonlinearity describing the system berge84 ; kodama91 . Different nonlinear power in the model imply different exponents, as described by Kodama et al. kodama91 . Kim et al. kim98  studied the exponent for the average length of laminar phases as function of the reinjection probability. Herein another average is explored to characterize the intermittency and its nonlinearity. It is the average of one variable of the system. Near bifurcation, approximate expressions for the natural measure, or probability density, are obtained for a map with type-III intermittency, in section II. Critical exponents are defined from analytical approximations and shown to be useful in characterizing the bifurcation.

Numerical experiments with the maps are presented in section III, verifying the proposed exponents. To test in a real physical system, a nonlinear circuit with a diode, similar to the one used in the early demonstrations of chaotic universal properties jeffries82 ; linsay81 ; testa82 ; buskirk85 , was set up as described in section IV and used to verify the exponents. The voltage across the diode, which is a dynamical variable of the system, was simultaneously measured in time series and average. The type-III intermittency bifurcation is well characterized both, using the experimental next peak value map and the average.

## II Normal Form map with type-III intermittency

To establish the new critical exponents for the averages one begins with the normal form of the map that has type-III intermittency berge84

$$M(X)=(1+ϵ^{})X+\alpha X^2+\eta ^{}X^3$$
(1)
The bifurcation, when $`X=0`$ ceases to be a stable fixed point, occurs at $`ϵ^{}=0`$. The second application of this map, in the approximation of small $`ϵ^{}`$ and X is given by

$$M^2(X)=(1+2ϵ^{})X+\eta X^3$$
(2)
where $`\eta =2(\eta ^{}+\alpha ^2)`$. If this coefficient is positive and $`ϵ^{}>0`$ one has type-III intermittency around $`X=0`$, provided a reinjection mechanism is introduced in eqs. (1) and (2). When $`\alpha 0`$ the map of eq. (1) is nonsymmetrical in $`X`$. Generally, any odd exponent in the nonlinear term of eq. (2) leads to the intermittency. Thus, a map (second iterate) defined in the interval $`[0,1]`$ with the type-III intermittency is manneville90

$$X_{n+1}=\left[(1+ϵ)X_n+X_n^z\right]\mathrm{mod}1$$
(3)
The value of $`z3`$ describes the nonlinear dependence of the subcritical bifurcation at $`ϵ=0`$. For $`ϵ1`$ the fixed point $`X=0`$ is unstable, but many iterates of the map, shown in Figure 1, falls near zero. These are called the laminar phase iterates.

For small $`ϵ`$ the natural measure is obtained, following the steps of Manneville manneville90a , as

$$\mu (ϵ,z,X,dX)=\left[_0^1\frac{dX}{ϵ+X^{(z1)}}\right]^1\frac{dX}{ϵ+X^{(z1)}}$$
(4)
Using this measure, the average of $`X`$ is

$$X(ϵ,z)=\left[_0^1\frac{dX}{ϵ+X^{(z1)}}\right]^1_0^1\frac{XdX}{ϵ+X^{(z1)}}$$
(5)
Its dependence on $`ϵ`$ and $`z`$ can be analytically established for small $`ϵ`$.

The simple case $`z=3`$ reduces to

$$X(ϵ,z=3)\frac{1}{\pi }ϵ^{1/2}\mathrm{ln}ϵ$$
(6)
For $`z>3`$ the general asymptotic expression is

$$X(ϵ,z)ϵ^\nu $$
(7)
where $`\nu `$ is the new exponent, whose value is

$$\nu =1/(z1)$$
(8)
This is the exponent that can be extracted from numerical and experimental systems to obtain the value of $`z`$. Similar expressions can be obtained for the second iterated map with negative values of $`X`$. If the original first iterated map, giving eq. (3) is fully odd in $`X`$ the total average must be zero. Once the reinjection or the map is not symmetrical, i. e., the value of some even nonlinear power (less than $`z`$) is non zero, the exponents of the averages can be obtained from the second iterate maps as eq. (3). A detailed study of these symmetry properties will be published elsewhere hugo2001 . One must recall that, in type-III intermittent chaos, an exponent already exists in the literature for the average length of laminar phases, given by kodama91 ; ono95  $`\beta =(z2)/(z1)`$. The simple relation exists, giving $`\nu +\beta =1`$. Therefore, as $`z`$ varies through $`3,5,\mathrm{}`$ the relative variation of these exponents are $`\mathrm{\Delta }\nu /\nu >\mathrm{\Delta }\beta /\beta `$. This means equal or better sensitivity in the determination of $`\nu `$ with respect to $`\beta `$. The most important fact is that those two exponents should be searched for independently in any experiment.

## III Numerical Averages and exponents in the map

To test the predictions of eq. (8) numerically the map of eq. (1) was restricted to $`|X|<1`$ by a modulo one reinjection: When the iterate gives $`X>1`$ or $`X<1`$ the integer part is subtracted. Averages were obtained directly from iterates of the map starting from random initial condition. A transient of $`2\times 10^4`$ iterates was eliminated from a total of $`10^6`$ calculated values at each of $`10^3`$ steps of $`ϵ^{}`$, between zero and $`10^3`$. For the average to be non null the value of $`\alpha `$ has to be different of zero. This breaks the symmetry between positive and negative values of X. Once the average is non null, its behavior for small $`ϵ`$ follows eq. (6). This is shown in Figure 2. The $`\alpha =0`$ case, where the average is always zero, is not represented.

The same slope is obtained for the averages with different values of $`\alpha `$ and $`\eta ^{}`$. The resulting averages have the same dependence in $`ϵ`$ as the one obtained from the map of second iterates (eq. (3)), calculated with the same number of iterates; all coinciding with the expression of eq. (6).

Figure 3 shows numerically calculated averages obtained directly from iterates of the map eq. (3), when $`z=3`$ and $`z=7`$. Again the initial condition at each value of $`ϵ`$ was taken at random. A transient of $`10^5`$ iterates was eliminated from a total of $`2\times 10^7`$ calculated values at each of $`10^3`$ steps of $`ϵ`$, between zero and $`10^1`$.

For $`ϵ<10^2`$ the $`z=3`$ and $`7`$ results are very well superimposed by the analytical curves made from equation (6), with $`\nu =0.50`$, and equation (7), with $`\nu =1/60.17`$, respectively. The two behaviors are clearly distinguished in the log-log plot. For the purpose of comparison with the already established exponents one should notice that the exponents for the average laminar phase are given in kodama91  as $`(z2)/(z1)`$ that is $`\beta =0.50`$ and $`\beta =0.83`$, respectively, verifying $`\nu +\beta =1`$.

Up to this point all results have concerned discrete maps. Dynamical fluxes with associated maps having intermittent chaos show similar critical exponents for the average of its continuous variables. This has been verified numerically with tangent bifurcation of type-I intermittency hugomapaave2000 ; hugoave2000  and will be demonstrated here for a physical chaotic oscillator. For the type-III intermittency the map extracted from a flux has the form of eq. (3) for its second iterates. A discretely sampled continuous flux may be viewed as a finite set of discrete stroboscopic maps, all with the same stroboscopic frequency but stroboscopic phases ranging from 0 to $`2\pi `$ in uniformly separated steps. The number of such maps gives the ratio of sampling frequency to stroboscopic frequency. The time average of the continuous flux is the average of the time averages of these maps. If a map from this series is nonsymmetrical, only in very special cases it would happen that this asymmetry be exactly canceled by the summation on the other maps. Thus, in general, its expected that if the flux leads to a nonsymmetrical map (no odd symmetry with respect to the unstable fixed point) the average of its continuous variable will exhibit the same exponent $`\nu `$ of eq. (8).

## IV Experimental exponents in a nonlinear circuit with intermittency

The physical experimental system to test the above exponents consisted of an RLC series circuit, where a nonlinear capacitance was implemented by a p-n junction diode testa82 ; buskirk85 . The circuit, presented in figure 4, was driven by a frequency and amplitude controllable oscillator with output impedance of $`50`$ Ohms. An inductor of $`0.1`$ Henry, an external resistor of $`13`$ Ohms and a (typical 1N4007) diode formed the circuit, whose linear oscillation regime had resonance frequency at $`150`$ kHz.

Dynamical bifurcations were produced by scanning the external drive frequency. Time series of the value of the voltage across the diode and the current in the circuit were collect with a $`12`$ bits resolution converter. The sampling rate was $`10^7`$ sample/s. Thus $`200`$ points were saved on each oscillator cycle. Capturing $`10^4`$ points at each value of the control frequency, a simple software searched for the maxima in the series. In wide range scans the typical bifurcation diagrams given by the peak value of the voltage across the diode, show clearly the well known results of period doubling cascades, chaotic windows and tangent bifurcations from chaos into periodic windows. jeffries82 ; testa82 ; buskirk85 ; ono95 .

A specific bifurcation, with intermittency and no bistability, was found and studied, scanning the external oscillator frequency from $`48708`$ Hz to $`48768`$ Hz by $`300`$ equal small steps. The drive voltage amplitude was fixed at $`2.7`$ V. The circuit was mounted inside an isolating box to prevent thermal drift effects donoso97 . The value of a control parameter $`ϵ`$ is obtained as the frequency detuning step divided by the critical frequency of the bifurcation, $`48762`$ Hz. Thus $`ϵ`$ varied in the range of $`10^4`$$`10^3`$. A segment of the voltage pulses is shown in figure 5. The signature of type-III intermittency is observed with the laminar events corresponding to the uniform oscillations alternating their peak value around 1V.

Multibranch maps constructed from the maxima of the voltage across the diode, extracted from series with $`6\times 10^5`$, could be approximated to one-dimensional ones, expected in the limit of infinite dissipation. Second return maps and laminar phase histograms, not shown here, indicate type-III intermittency. Figure 6 shows the average length of laminar phases for different values of $`ϵ`$. Each histogram used was obtained with $`6\times 10^5`$ points. A theoretical fitting kodama91  is best with an exponent $`\beta =0.62`$. Notice that for the map the predictions are $`\beta =0.5`$ for $`z=3`$ and $`\beta =0.75`$ for $`z=5`$. However the experimental data always has excess of laminar events identified with short length griffith97  and this effect gives a bigger experimental value for $`\beta `$ in the fittings. Therefore, $`z=3`$ is the best odd value.The excess of short laminar phase events may be related to the finite dissipation rate in the phase space of the experimental system and the consequent non unidimensional maps and also to a nonuniform density of probability for the reinjection, as discussed by Kim et al. kim98 .

The peak value of the voltage, which gave the second return map and the histograms is shown in Figure 7(a). The laminar phases are the oscillations with repetitive visits to the maximum value near $`1V`$ in the figure. As the segments represented in figure 5(a) are short in time for many values of parameters they show a single laminar event with all points accumulated around the unstable orbit. Also shown is the simultaneously acquired average of the voltage across the diode. It was obtained with a simple electronic integrator, having a time constant of $`3`$ seconds. To account for long laminar phase events and (which is equivalent) decrease the average fluctuations the scan lasted 50 minutes.

The experimental average voltage was fitted to the expression

$$X(ϵ)=ϵ^\nu \mathrm{ln}ϵ$$
(9)
The exponent $`\nu =0.55`$ gave an excellent agreement with the experimental plot, as shown in figure 7. Attempts of fittings to higher values of $`z`$, i.e., to eq. (7), failed. It is worth noticing that the value predicted for $`z=5`$ is $`\nu =0.25`$. Thus, the experimental average consistently verifies $`z=3`$ for the nonlinearity of this bifurcation in the circuit. The confidence for the experimental values of $`\beta `$ and $`\nu `$ from fittings using standard $`\chi ^2`$ procedure is better than 2%. Therefore, the result $`\beta +\nu =1.17\pm 0.02`$ shows that a discrepancy remains between the experiments and the unidimensional map model. While the bigger value obtained for $`\beta `$ has been attributed to an excess of short laminar phase events kim98  no such study of deviations exists for the exponent in the averages.

## V Conclusion

In conclusion, critical exponents for the averages of one dynamical variable are established analytically for type-III intermittent chaotic maps. Those exponents are directly related to the nonlinear power law of the normal form of the maps. They have a simple relation with the exponents of the average length of laminar iterates in the same systems. All these properties were verified in numerical experiments with maps.

Physical experiments with a continuous flux were also done to demonstrate the exponents in averages. The average voltage across a diode in a chaotic electronic circuit was measured while the drive frequency was scanned through bifurcations. A bifurcation from chaos into periodic pulsation, shown to be type-III intermittent, gives an exponent in agreement with the cubic nonlinearity. This result is consistently verified in the exponents of the average length of laminar phases, extracted from histograms of the peak pulse voltages in the circuit.

Averages of dynamical variables have been proposed to get the signature of the Lorenz chaos bifurcation lawandy87 , have been numerically studied in critical bifurcations rajasekar00 , and experimentally measured in bifurcating pulsed lasers oliveira96 . However no systematic study of critical exponents have been done.

The technique of measuring averages of dynamical variable in intermittent chaos is a complementary procedure to investigate bifurcations of nonlinear systems. The experimental average may also be advantageous when detection noise for an specific variable has a bandwidth overlapping the frequency bandwidth of the chaotic oscillations. For systems with high frequency noise it is naturally bound to be more sensitive. Its experimental motivation is therefore enhanced to characterize ultrafast chaotic oscillators, like diode lasers, using slow time detection techniques. One extension relevant to the work presented here is the study of the exponents for the averages in bifurcations with nonuniform reinjection in the intermittency, as studied by Kim et al. kim98 . A relation between exponents of the averages and the exponents of the average length of laminar phases events should exist generalizing the $`\nu +\beta =1`$ given here. Another potential use for the averages near bifurcation would be to complement the confrontation between experimental data and model as extracted from nonlinear data analysis hegger98 ; timmer2000 . The models inferred from the data analysis must have bifurcations consistent with the experimental system. These bifurcations could be tested comparing both the exponents of laminar events and the exponents of averages of dynamical variables. The Critical exponents have been introduced for many types of bifurcations in chaotic systems grebogi86 . Their presence in simple, experimentally accessible, statistical properties, as the averages and its higher moments, are under investigation. The earliest citation of averages in chaos can be traced to the original propositions of unpredictability in deterministic chaos lorenz64 .

## acknowledgments

Work partially supported by Brazilian Agencies: Conselho Nacional de Pesquisa e Desenvolvimento (CNPq) and Financiadora de Estudos e Projetos (FINEP)


text:
# UNUSUAL SYMMETRIES IN THE KUGEL-KHOMSKII HAMILTONIAN

## Abstract

The Kugel-Khomskii Hamiltonian for cubic titanates describes spin and orbital superexchange interactions between $`d^1`$ ions having three-fold degenerate $`t_{2g}`$ orbitals. Since orbitals do not couple along “inactive” axes, perpendicular to the orbital planes, the total number of electrons in $`|\alpha `$ orbitals in any such plane and the corresponding total spin are both conserved. A Mermin-Wagner construction shows that there is no long-range spin ordering at nonzero temperatures. Inclusion of spin-orbit coupling allows such ordering, but even then the excitation spectrum is gapless due to a continuous symmetry. Thus, the observed order and gap require more symmetry breaking terms.

High temperature superconductivityB  and colossal magnetoresistanceCMR  sparked much recent interest in the magnetic properties of transition metal oxides, particularly those with orbital degeneracyKKrev ; TN . In many transition metal oxides, the $`d`$ electrons are localized due to the large on-site Coulomb interaction, $`U`$. Assuming a simple Hubbard model, with a typical nearest neighbor (nn) hopping energy $`t`$, the low energy behavior can be described by an effective superexchange model, which involves only nn spin and orbital coupling, with energies of order $`ϵ=t^2/U`$. In cubic oxide perovskites, the crystal field of the surrounding oxygen octahedra splits the $`d`$-orbitals into a two-fold degenerate $`e_g`$ and a three-fold degenerate $`t_{2g}`$ manifold. In most cases, these degeneracies are further lifted by a cooperative Jahn-Teller (JT) distortionKKrev , and the low energy physics is well described by an effective superexchange spin-only model ANDERSON ; TY1 ; TY2 . However, some cubic perovskites, such as the titanates (RTiO3, where R= La, Y, etc), have only a small JT distortionJTdist , in spite of the orbital degeneracynoJT . This distortion, as well as the small orthorhombic deviation from the cubic symmetry, were not even observed in Refs. LTO, ; YTO, . Since these distortions are small, several theoretical papers chose to ignore them and assume cubic symmetry. The corresponding cubic model has been taken as the “minimal” model needed to explain the physics in these materials. In the present Letter we show that although this model is of great theoretical interest, it is insufficient to explain the experiments.

For the cubic titanates, there is one $`d`$ electron in the $`t_{2g}`$ degenerate manifold, which contains the wavefunctions $`|Xd_{yz}`$, $`|Yd_{xz}`$, and $`|Zd_{xy}`$. The large degeneracy of the resulting ground states, which involve both the spin and the orbital degrees of freedomKKrev ; TN ; KK , may then yield rich phase diagrams, with exotic types of order, involving a strong interplay between the spin and orbital sectors (e.g. TN, ; LTO, ; YTO, ), justifying the broad theoretical interest in this cubic limit. As we show, the corresponding superexchange Hamiltonian (hereafter called the cubic Kugel-Khomskii (KK) modelKK ) contains several interesting hidden symmetries. In addition, our analysis shows that the KK Hamiltonian cannot yield some of the predictions which were claimed in the literature to follow from it. In particular, it has been suggested GK1  that the KK Hamiltonian gives rise to an ordered isotropic spin phase at non-zero temperatures, and that an energy gap in the spin excitations can be caused by spin-orbit interactions GK2 . We use the symmetries of the KK Hamiltonian to show that both of these predictions cannot hold. The observed long range order and finite gap LTO  must therefore be based on more complicated Hamiltonians, which go beyond the scope of this paper.

As can be seen from Fig. 1, cubic symmetry implies no hopping (via oxygen $`p`$-states) among orbitals of type $`|\alpha `$ along the $`\alpha `$-axis. Neglecting the direct Ti-Ti hopping, KK called this axis the “inactive axis” for $`\alpha `$ orbitals. This statement forms the basis for the remarkable symmetry properties of the KK Hamiltonian reported in this Letter. Apart from constant terms, the perturbative expansion of the Hubbard Hamiltonian with hopping of this type, to order $`ϵ`$, yields the cubic KK Hamiltonian, $`=_x+_y+_z`$, where

$`_\alpha =ϵ{\displaystyle \underset{ij\alpha }{}}{\displaystyle \underset{\beta ,\gamma \alpha }{}}{\displaystyle \underset{\sigma ,\eta }{}}c_{i,\beta ,\sigma }^{}c_{i,\gamma ,\eta }c_{j,\gamma ,\eta }^{}c_{j,\beta ,\sigma }.`$ (1)
and $`ij\alpha `$ denotes a nn bond along the $`\alpha `$-axis. Here, $`c_{i,\beta ,\sigma }^{}`$ creates an electron at site $`i`$ in a $`\beta `$ orbital with spin $`\sigma `$, and one assumes that there is exactly one electron on each site, i.e. $`_\alpha n_{i\alpha }=1`$, with $`n_{i\alpha }_\sigma c_{i,\alpha ,\sigma }^{}c_{i,\alpha ,\sigma }`$.

For some purposes, it is convenient to separate the spin and orbital degrees of freedomKK . Defining the spin of an electron at site $`i`$ as $`𝐒_i`$, one has

$`_\alpha =(ϵ/2){\displaystyle \underset{ij\alpha }{}}\left(1+4𝐒_i𝐒_j\right)J_{ij}^\alpha ,`$ (2)
$`J_{ij}^z=n_{ix}n_{jx}+n_{iy}n_{jy}+a_i^{}b_ib_j^{}a_j+b_i^{}a_ia_j^{}b_j.`$ (3)
Here, $`a_i^{}`$ and $`b_i^{}`$ create spin-less electrons in orbitals $`|X`$ and $`|Y`$, respectively, and $`n_{ix}=a_i^{}a_i`$ etc.

Both Eqs. (1) and (3) imply that whenever an $`\alpha `$-orbital electron is destroyed, an $`\alpha `$-orbital electron is created on either the same or another site. Therefore, the total number of electrons in each orbital is a good quantum number: any eigenfunction can be labeled by the total number of electrons in each orbital (i.e. $`N_X,N_Y,`$ and $`N_Z`$)IHM . Furthermore, when an $`\alpha `$-orbital electron is destroyed, it is replaced by another $`\alpha `$-orbital electron in the same plane perpendicular to the inactive ($`\alpha `$) axis. Thus, for the $`n`$’th plane perpendicular to the $`\alpha `$-axis, the total number $`N_{n\alpha }`$ of electrons in orbital $`|\alpha `$ is conserved, i. e. it is a good quantum number. For example, for a cube of eight sites, the numbers $`N_{1X}`$ and $`N_{2X}`$, which respectively are the numbers of X-orbital electrons in each of the two planes perpendicular to the $`x`$-axis, are conserved, and similarly for $`y`$ and $`z`$. Thus, the states of the cube can be labeled by the six quantum numbers, $`(N_{1X},N_{2X},N_{1Y},N_{2Y},N_{1Z},N_{2Z})`$.

Remarkably, there are more conserved quantities, associated with electron spins. Defining the spin of an electron in orbital $`\alpha `$ at site $`i`$ as $`𝐒_{i\alpha }_{\rho \eta }c_{i,\alpha ,\rho }^{}[\stackrel{}{𝝈}]_{\rho \eta }c_{i,\alpha ,\eta }/2`$, where $`\stackrel{}{𝝈}`$ represents the vector of Pauli matrices, and the total spin of all such electrons - located in an arbitrarily chosen plane #$`n`$ perpendicular to the inactive $`\alpha `$ axis - as $`\stackrel{}{𝒮}_{n\alpha }_{in}𝐒_{i\alpha }`$, we next perform a uniform, but arbitrary, rotation of all the $`𝐒_{i\alpha }`$’s with $`in`$: we introduce an arbitrary $`2\times 2`$ unitary matrix $`𝐔^{(𝐧)}`$, and write

$`c_{i,\alpha ,\sigma }^{}`$ $`=`$ $`{\displaystyle \underset{\eta }{}}U_{\sigma ,\eta }^{(n)}d_{i,\alpha ,\eta }^{},in.`$ (4)
Electrons in other orbitals or in other lattice planes are not affected by this transformation. Substitution of this transformation into Eq. (1) shows that it leaves $``$ invariant. As a consequence of this symmetry, if one assumes long-range spin order, the spins associated with $`\alpha `$-orbitals within any given plane can be rotated at zero cost in energy, thereby destroying the supposed long-range order. It also follows that $`\stackrel{}{𝒮}_{n\alpha }`$ commutes with $``$, and thus both $`|\stackrel{}{𝒮}_{n\alpha }|^2`$ and $`[\stackrel{}{𝒮}_{n\alpha }]_z`$ are good quantum numbers for each value of $`n`$ or $`\alpha `$. These symmetries can also be obtained from the original Hubbard model, provided one neglects Coulomb exchange interactions.

This situation allows a rigorous proof of the nonexistence of long-range spin order at any nonzero temperature for $``$ of Eq. (1). Following the procedure of Mermin and Wagner (MW)MERMIN , we choose

$`C`$ $`=`$ $`\widehat{S}_\alpha ^+(𝐤)={\displaystyle \underset{𝐑}{}}e^{i𝐤𝐑}c_{𝐑,\alpha ,}^{}c_{𝐑,\alpha ,},`$
$`A`$ $`=`$ $`\widehat{S}_\alpha ^{}(𝐤+𝐊)={\displaystyle \underset{𝐑}{}}e^{i(𝐤+𝐊)𝐑}c_{𝐑,\alpha ,}^{}c_{𝐑,\alpha ,},`$ (5)
where $`𝐊`$ is the wavevector of the order we wish to discuss. Here $`a𝐊=(\pi ,\pi ,\pi )`$ is most relevant. Assuming long range order of $`\widehat{S}_\alpha ^z(𝐊)s_{\alpha ,z}`$ and a corresponding staggered field $`h`$ (for spins in the $`\alpha `$-orbital), we end up with the MW-like bound

$`12kT|s_{\alpha ,z}|^2{\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤}{}}\left[hs_{\alpha ,z}+\widehat{J}_\alpha (𝐤)\right]^1,`$ (6)
where $`\widehat{J}_\alpha (𝐤)_{\stackrel{}{\delta }\alpha }(1e^{i𝐤\stackrel{}{\delta }})`$ is proportional to the $`𝐤`$-dependent parts in the Fourier transform of the non-zero nn spin-spin interaction (with nn vector $`\stackrel{}{\delta }`$) in the Hamiltonian. Since spins in orbital $`\alpha `$ couple only within planes perpendicular to the $`\alpha `$-axis, it follows that $`\widehat{J}_\alpha (𝐤)2_{\beta \alpha }[1\mathrm{cos}(k_\beta a)]_{\beta \alpha }a^2k_\beta ^2a^2k_{,\alpha }^2`$, with no dispersion in the $`\alpha `$-direction. For systems in $`d3`$ dimensions, the sum in Eq. (6) diverges as $`h0`$, implying that $`s_{\alpha ,z}`$ must go to zero. The conclusion is that the KK model is at its lower critical dimension $`d_<=3`$ and does not support long-range spin order at $`T>0`$. As we show elsewhere, a similar proof can be formulated for the original Hubbard modeldetails .

The same conclusion also follows from a renormalization group analysis of the model at finite $`T`$details . Generalizing to $`m`$ orbitals and $`n`$-component spins, the spin free energy functional maps onto that of the “canonical” $`nm`$-component spin problem DG , but with a $`(d1)`$-dimensional transverse gradient term,

$`F`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐪\alpha }{}}(r+q_{,\alpha }^2)\widehat{𝐒}_\alpha (𝐪)\widehat{𝐒}_\alpha (𝐪)`$ (7)
$`+`$ $`{\displaystyle \underset{𝐑}{}}\left(u{\displaystyle \underset{\alpha }{}}|\widehat{𝐒}_\alpha (𝐑)|^4+v{\displaystyle \underset{\alpha <\beta }{}}|\widehat{𝐒}_\alpha (𝐑)|^2|\widehat{𝐒}_\beta (𝐑)|^2\right),`$
where $`\widehat{𝐒}_\alpha (𝐪)`$ is the Fourier transform of $`𝐒_\alpha (𝐑)`$. Similar forms arise in connection with Lifshitz-like behaviorAARG . This anisotropic gradient term shifts both the upper and the lower critical dimensions up by 1. For $`3<d<5`$ dimensions and $`n>1`$ it yields decoupled $`n`$-component critical behavior. The free energy (7) also reflects the symmetry with respect to independent rotations of the spin $`𝐒_\alpha `$ associated with the single orbital $`\alpha `$.

These symmetries and conservation laws are very useful in the exact numerical diagonalization of finite Ti clusters, which indeed confirms their validity. We demonstrate this for a cube of 8 sites (Fig. 2). Since the Hamiltonian commutes with the total spin $`𝐒=_i𝐒_i`$, the eigenstates can be identified by the quantum numbers $`S`$ \[where $`𝐒^2=S(S+1)`$\] and $`S_z`$. Since the energy does not depend on $`S_z`$, it suffices to study the subspace of $`70\times 3^8=459270`$ states with $`S_z=0`$. A numerical analysis of the low-energy spectrum of this huge sparse matrix yielded 3 degenerate $`S=0`$ ground states, $`\mathrm{\Psi }_x`$, $`\mathrm{\Psi }_y`$, and $`\mathrm{\Psi }_z`$, related by cyclic permutations. $`\mathrm{\Psi }_z`$ has $`N_x=N_y=4`$. This information suffices to find the manifold containing each of these ground states. Since $`\mathrm{\Psi }_z`$, for instance, is not degenerate within the manifold $`N_x=N_y=4`$, it must also have the quantum numbers $`N_{1X}=N_{2X}=2`$, $`N_{1Y}=N_{2Y}=2`$, and $`N_{1Z}=N_{2Z}=0`$. A nonsymmetric choice like $`N_{1X}=3`$ and $`N_{2X}=1`$ would be degenerate with $`N_{1X}=1,N_{2X}=3`$. The lack of degeneracy also implies that the total spin of the $`N_{n,\alpha }=2`$ electrons in orbital $`\alpha `$ in the $`n`$’th plane perpendicular to the $`\alpha `$-axis must be $`𝒮_{n\alpha }=0`$. Examples of such configurations, containing dimers of $`\alpha `$-orbital electrons in the $`\alpha `$-planes, are shown in Fig. 2. The Hamiltonian allows an exchange of an $`X`$-electron with a $`Y`$-electron only along the $`z`$-axis (the only axis along which both types can hop). Starting from Fig. 2(a), and performing all such possible exchanges, creates a manifold of 16 states (3 of which are shown in the figure). Two other states with the same $`N_{n\alpha }`$’s, but with the dimers along the $`z`$-axis, form another manifold, of higher energy. Indeed, a diagonalization of the resulting $`16\times 16`$ matrix reproduced the same ground state energy as found from the $`459270\times 459270`$ matrix, demonstrating the power and the correctness of these symmetries. We are currently extending these numerical studies to even larger systems (such as $`N=16`$ sites) to better understand the nature of the ground state in a real system.

Since the KK Hamiltonian (1) forbids long range spin order at $`T>0`$, the existence of such order (as in LaTiO3 LTO ) must result from some additional mechanism. Even for cubic symmetry, such mechanisms could include the small direct Ti-Ti hopping along the inactive axis, Coulomb exchange terms in the original Hubbard model, or the spin-orbit interaction. In the real orthorhombic titanates one must also include JT distortions and oxygen octahedra rotations. A full discussion of all these effects lies beyond the present paper. Since the present paper concerns mainly symmetry arguments, we concentrate here on adding the spin-orbit interaction to the cubic KK Hamiltonian, where we can use such arguments to show the absence of a spin gap. Specifically, the spin-orbit Hamiltonian is

$`_{SO}=\lambda {\displaystyle \underset{i\alpha \beta \sigma \sigma ^{}\mu }{}}L_{\alpha \beta }^\mu c_{i,\alpha ,\sigma }^{}[𝝈_\mu ]_{\sigma \sigma ^{}}c_{i,\beta ,\sigma ^{}},`$ (8)
where $`L_{\alpha \beta }^\mu \alpha |L^\mu |\beta `$ is the orbital angular momentum matrix element. Since $`_{SO}`$ mixes orbitals, an $`\alpha `$-electron can hop via $`_{SO}`$ to orbital $`\beta `$, then hop to a $`\beta `$-orbital on a nn along the $`\alpha `$-axis, and finally use $`_{SO}`$ to return to orbital $`\alpha `$. This generates an effective hopping between $`\alpha `$-orbitals along the “inactive” $`\alpha `$-axis, invalidating the above arguments, shifting the lower critical dimension for total spin ordering back to $`d_<=2`$ and restoring long range spin order at $`d=3`$. This mixing also eliminates the independent symmetries which we found for electrons within each orbital separately. However, as discussed below, there still remain some global symmetries for the total spin. Based on the signs of the leading couplings, we assume that the total spin orders antiferromagnetically, and proceed to show that the spin wave excitations in the ordered phase must be gapless.

Again, an exact symmetry analysis clarifies the situation. For electrons within each of the three degenerate $`t_{2g}`$ orbitals discussed here, we introduce the following canonical transformationTY1  from spin to pseudo-spin:

$`c_{i,\alpha ,\sigma }^{}`$ $`=`$ $`{\displaystyle \underset{\eta }{}}𝐔_{\sigma ,\eta }^{(\alpha )}d_{i,\alpha ,\eta }^{},`$ (9)
where $`U^\alpha =𝝈_\alpha `$ represents a different rotation for spins in different orbitals. As discussed in Ref. TY1, , all the terms in both the original $``$ and in $`_{SO}`$ now contain only combined operators of the form $`_\mu d_{i,\alpha ,\mu }^{}d_{i,\beta ,\mu }`$, with coefficients which do not depend on the pseudo-spin indices $`\mu `$ (see Eq. (6) in Ref. TY1, ). These combined operators, and therefore also the full Hamiltonian, are rotationally invariant in pseudo-spin space. Said differently, the Hamiltonian is invariant with respect to a transformation on the original spins of the form $`\stackrel{~}{c}_{i,\alpha ,\mu }^{}=_\mu ^{}V_{\mu ,\mu ^{}}^{(\alpha )}c_{i,\alpha ,\mu ^{}}^{}`$ with $`𝐕^{(\alpha )}=𝝈_\alpha 𝐔𝝈_\alpha `$, where $`𝐔`$ is an arbitrary unitary matrix. Thus the system possesses a continuous symmetry, but it is not the usual symmetry with respect to rotation of the total spin. In the antiferromagnetically ordered phase, the spin staggered moment selects an orientation, and therefore the pseudo-spin will also exhibit broken symmetry. Rotation of the pseudo-spin will give rise to a manifold of zero energy states. This continuous symmetry guarantees that we have a (probably propagating) zero-energy hydrodynamic modeTCL . The rigorous conclusion is then that spin-orbit interaction permits the existence of long-range order at nonzero temperatures, but does not cause a gap in the elementary excitation spectrum, contrary to the assertion in Ref. GK2, . Since our argument is based on symmetry considerations, it holds no matter what type of fluctuation is considered, and regardless of the orbital ordering (long ranged or liquid). In analogy with results of Refs. TY1, ; TY2, ; BS, , it is probable that when Coulomb exchange interactions and/or canting of the Ti–O–Ti bonds are included, spin-orbit interactions would lead to an energy gap in the excitation spectrum.

In conclusion, we uncovered several novel symmetries of the KK Hamiltonian for cubic $`t_{2g}`$ systems. It is surprising that the KK Hamiltonian has been widely used in the study of interesting spin-orbital physics of transition metal oxides for a long time but yet its remarkable symmetry properties were missed until now. Using these symmetries, we rigorously showed that the KK Hamiltonian without spin-orbit interactions does not permit the development of long-range spin order in a three dimensional cubic lattice at nonzero temperature. Inclusion of spin-orbit interactions allows the formation of long-range spin order, but the excitation spectrum is gapless.

ABH thanks NIST for its hospitality during several visits when this work was done. We acknowledge partial support from the US-Israel Binational Science Foundation (BSF). The TAU group is also supported by the German-Israeli Foundation (GIF).