text:
# Plasma formation from ultracold Rydberg gases

## Abstract

Recent experiments have demonstrated the spontaneous evolution of a gas of ultracold Rydberg atoms into an expanding ultracold plasma, as well as the reverse process of plasma recombination into highly excited atomic states. Treating the evolution of the plasma on the basis of kinetic equations, while ionization/excitation and recombination are incorporated using rate equations, we have investigated theoretically the Rydberg-to-plasma transition. Including the influence of spatial correlations on the plasma dynamics in an approximate way we find that ionic correlations change the results only quantitatively but not qualitatively.

Advances in cooling and trapping of neutral gases have opened up a new branch of atomic physics, namely dynamics in ultracold ($`T1`$ K) systems. One interesting topic is the physics of ultracold Rydberg gases and ultracold neutral plasmas. In recent experiments, ultracold neutral plasmas have been produced from a small cloud of laser-cooled atoms confined in a magneto-optical trap Killian et al. (1999); Kulin et al. (2000); Killian et al. (2001); Robinson et al. (2000); Eyler et al. (2000). In the experiments performed at NIST Killian et al. (1999); Kulin et al. (2000); Killian et al. (2001), a plasma was produced by photoionizing laser-cooled Xe atoms with an initial ion temperature of about $`10\mu \mathrm{K}`$. By tuning the frequency of the ionizing laser, the initial electron energy $`E_e`$ was varied corresponding to a temperature range $`1\mathrm{K}<E_e/k_\mathrm{B}<1000\mathrm{K}`$, and the subsequent expansion of the plasma into the surrounding vacuum was studied systematically. Remarkably, a significant amount of recombination was observed as the plasma expands, leading to the formation of Rydberg atoms from the plasma. In a complementary type of experiment Robinson et al. (2000); Eyler et al. (2000), ultracold atoms were laser-excited into high Rydberg states rather than directly ionized. In these experiments, also the reverse process has been observed, namely the spontaneous evolution of the Rydberg gas into a plasma.

An issue in the theoretical considerations Kuzmin and O’Neil (2002); Robicheaux and Hanson (2002); Mazevet et al. (2002); Tkachev and Yakovlenko (2001) stimulated by the experiments has been the question whether the evolving plasma is strongly coupled or not. Due to the very low initial temperature, the electron Coulomb coupling parameter $`\mathrm{\Gamma }_e(t=0)`$ is found to be significantly larger than one ($`\mathrm{\Gamma }_e=e^2/(ak_\mathrm{B}T_e)`$, where $`a`$ is the Wigner-Seitz radius). On the other hand, it has been pointed out in Kuzmin and O’Neil (2002) that, for the initial conditions of the NIST experiments Killian et al. (1999); Kulin et al. (2000); Killian et al. (2001), the development of equilibrium electron-electron correlations leads to a rapid heating of the electron gas, which brings $`\mathrm{\Gamma }_e`$ down to order unity on a very short timescale. In Kuzmin and O’Neil (2002), this was demonstrated by a molecular-dynamics simulation of the electron and ion motion. The calculation was limited to a short time interval ($``$ 100 nanoseconds) in the initial stage of the plasma evolution due to the large numerical effort required. For the first quantitative comparison with experiment the plasma dynamics has been modeled within a kinetic approach, while ionization, excitation and recombination has been treated by a separate set of rate equations Robicheaux and Hanson (2002). The electron dynamics was described in an adiabatic approximation with the atom and ion temperatures set to zero. Since this model does only account for the mean-field potential created by the charges, possible correlation heating could not be described. Nevertheless, due to heating by three-body recombination events, $`\mathrm{\Gamma }_e`$ does not exceed a value of $`0.2`$ during the plasma expansion. Hence, the influence of electronic correlations on the plasma dynamics could be neglected. With zero ionic temperature the ionic Coulomb coupling parameter, however, is infinite in the framework of this model. Consequently, the role of ion-ion correlations could not be estimated.

Our description is similar to that of Robicheaux and Hanson (2002) with small differences, e.g., inclusion of black-body radiation as a source for photoionization. This is necessary to describe the initial ionization of the Rydberg gas in the experiments Robinson et al. (2000); Eyler et al. (2000) since the ionization rate through cold atom-atom collisions is much lower at these densities. More importantly, in addition to Robicheaux and Hanson (2002) we allow for possible spatial correlation effects. Briefly, we use a set of kinetic equations for the plasma expansion, combined with rate equations for the description of ionization/excitation and recombination. The kinetic equations are derived from the first equation of the BBGKY hierarchy, which yields the standard Vlasov equation for the evolution of a collisionless plasma, augmented by additional terms accounting for spatial correlations. Since the electronic Coulomb coupling parameter is well below unity under all experimental conditions realized so far (as we have convinced ourselves including electronic correlations on the basis of two-component Debye-Hückel theory; see also Robicheaux and Hanson (2002)) we neglect them and treat the ions as a one-component plasma embedded in a neutralizing background.

Assuming that the correlation function vanishes at distances larger than the correlation length $`a_c`$ and that the spatial density varies slowly over a distance $`a_c`$, we find that the force induced by ionic correlations can be expressed in terms of the local correlation energy per particle

$$u_{ii}(𝐫)=\frac{e^2}{2}\rho _i(𝐫)𝑑𝐲\frac{g_{ii}(y;𝐫)}{y},$$
(1)
where the spatial density $`\rho _i`$ is related to the one-particle distribution function $`f_i`$ by $`\rho _i=𝑑𝐯f_i`$ and $`g(y)`$ is the correlation function. The correlation energy for a one-component plasma in thermal equilibrium is a function of the Coulomb coupling parameter only Dubin and O’Neil (1999), and simple analytical formulae have been given for a wide range of $`\mathrm{\Gamma }`$. We use

$$u_{ii}^{(\mathrm{eq})}=\frac{9}{10}k_\mathrm{B}T_i\mathrm{\Gamma }_i=\frac{9e^2}{10}\left(\frac{4}{3}\pi \rho _i\right)^{1/3},$$
(2)
derived from the ion sphere model which is a good approximation for $`\mathrm{\Gamma }>1`$ Ichimaru (1982).

The timescale for relaxation towards the individual equilibrium of electrons and ions, respectively, is determined by the corresponding inverse plasma frequency $`\omega _\mathrm{p}^1=\sqrt{m/4\pi \rho e^2}`$ Murillo (2001). Under typical conditions of the experiments Killian et al. (1999); Kulin et al. (2000); Robinson et al. (2000); Killian et al. (2001), it is of the order of $`10^9`$s for electrons, while equilibration of the ions takes between $`10^7`$s and $`10^5`$s depending on $`\rho `$. These timescales for equilibration have to be seen relative to the expansion time of the plasma, which is of the order of several $`10^6`$s. Hence, an adiabatic approximation for the electron distribution function can be safely applied, such that the kinetic equation for the electronic component of the plasma yields a relation allowing one to express the mean-field potential in terms of the density $`\rho _e`$, which can be inserted into the corresponding equation for the ions.

The much slower equilibration of the ions renders an adiabatic approximation a priori difficult. However, even in situations where static parameters, e.g., relevant masses, periods or rates, speak against an adiabatic treatment it is sometimes made possible through dynamical adiabaticity (see, e.g., Tanner et al. (2000)). Of course, this can only be verified a posteriori by more elaborate calculations without adiabatic approximation. For now, we simply assume an adiabatic time evolution of the plasma, i.e. we estimate the ionic correlation energy from its equilibrium value eq. (2) and use a homogeneous ion temperature at all times. This is the opposite limit to a zero ionic temperature which would never lead to equilibration. Hence, our results are expected to clarify whether (and to what extent) correlation effects can have any influence on the dynamics at all.

The large difference in timescales between ionic and electronic motion justifies the quasineutral approximation $`\rho _e=\rho _i=\rho `$ Dorozhkina and Semenov (1998) applicable under the experimental conditions in Killian et al. (1999); Kulin et al. (2000); Robinson et al. (2000); Killian et al. (2001). Under these circumstances a closed equation for the ion distribution function is obtained and the resulting equation without the correlation term permits selfsimilar analytical solutions Dorozhkina and Semenov (1998). One of them is the gaussian profile describing the initial state in the experiments under consideration. Exact selfsimilar solutions exist even in the case of an additional linear force Dorozhkina and Semenov (1999). This condition is satisfied by the correlation force to a good approximation save for the outer periphery of the plasma Pohl et al. . Hence, we can write

$$f_i\mathrm{exp}\left(\frac{r^2}{2\sigma ^2}\frac{m_i\left(𝐯𝐰(𝐫)\right)^2}{2k_\mathrm{B}T_i}\right),$$
(3)
where $`\sigma `$ is the width of the spatial density distribution and $`𝐰(𝐫)=\gamma 𝐫`$ is the hydrodynamic velocity (see alsoRobicheaux and Hanson (2002)). Substitution of Eq. (3) finally leads to the set of equations

$`{\displaystyle \frac{d\sigma ^2}{dt}}`$ $`=`$ $`2\gamma \sigma ^2`$ (4)
$`{\displaystyle \frac{d\gamma }{dt}}+\gamma ^2`$ $`=`$ $`{\displaystyle \frac{N_e}{M\sigma ^2}}\left[k_\mathrm{B}\left(T_e+T_i\right)+W_c\right],`$
where $`W_c=1/3𝑑𝐫\rho (u_{ii}\rho )/\rho `$ arises from the correlation pressure and $`M`$ is the total mass of the plasma. The thermal energy is determined from the total energy of the system

$$E_{\mathrm{tot}}=\frac{3}{2}N_ek_\mathrm{B}\left(T_e+T_i\right)+\frac{3}{2}M\gamma ^2\sigma ^2+U_c,$$
(5)
i.e. the sum of the kinetic energy and the correlation energy $`U_c=𝑑𝐫u_{ii}\rho `$. The set of equations (4) reduces to the one used in Robicheaux and Hanson (2002) if the correlation term is dropped and $`T_i0`$. More recently, it has been shown by comparison with more sophisticated calculations that this simple ansatz provides a good description of the plasma dynamics even if some of the assumptions, such as the quasineutrality of the plasma and the gaussian distribution of the densities, only hold approximately Robicheaux and Hanson (2003) .

The low-temperature enhancement of the expansion velocity of the ultracold plasma produced in Kulin et al. (2000) can be well described by combining a hydrodynamic description of the expansion dynamics with conventional rate equations accounting for inelastic collisions between the plasma particles and Rydberg atoms Robicheaux and Hanson (2002). As a result, the total kinetic energy becomes a function of time since the Rydberg atoms act as energy sinks and sources. In our calculations, we include bound-bound transitions by electron impact excitation and deexcitation, and bound-free transitions by three-body recombination, electron impact ionization, and black-body radiation. We use the collision rate coefficients derived by Mansbach and Keck Mansbach and Keck (1969) and black-body photoionization is described by first-order perturbation theory Gallagher (1994), using an asymptotic expression for the atomic oscillator strengths Johnson (1972).

Finally, the question of the expansion velocity of the Rydberg atoms has to be addressed, since they are not driven by the Coulomb drag of the expanding electrons. Under typical experimental conditions ($`\rho _e=10^9\mathrm{cm}^3`$, $`T_e=20`$ K), the timescale for electron-atom-collisions is of the order of $`10^8`$s. Therefore, in the course of the expansion of the plasma ($`10^5`$s) any given atom with binding energy of the order of several $`k_\mathrm{B}T`$ will recombine and re-ionize many times, changing its character between ion and Rydberg atom. In addition, collisions between ions and atoms significantly equilibrate the hydrodynamical velocities even for lower Rydberg states Nemirovsky et al. (2002). Following this reasoning we assume equal hydrodynamical velocities and density profiles for the ions and atoms. This implies that the expansion of the neutral Rydberg atoms can be simply taken into account by changing the mass of the ions to an effective mass $`M`$ in Eq. (4) which is the total mass of the total system (plasma + atoms).

We have simulated the expansion of an initially fully ionized plasma as reported in Kulin et al. (2000); Killian et al. (2001). In accordance with previous calculations Robicheaux and Hanson (2002), we find quantitative agreement for the electron energy dependence of the asymptotic expansion velocity. Furthermore, the calculations qualitatively reproduce the nonmonotonic time dependence of the number of detected atoms observed in Killian et al. (2001). We note that the ionic correlations do not significantly change the plasma dynamics for these experiments.

The results of our calculations for the experiments Robinson et al. (2000); Eyler et al. (2000) are summarized in figures 1 to 3.

Figure 1 shows the degree of ionization of the Rydberg gas, i.e. the plasma fraction of the system, as a function of time during the expansion of the system. For comparison, the result without correlation ($`U_c=T_i0`$ in eqs. (4) and (5)) is also shown. As has been observed experimentally, it takes of the order of two to three microseconds before significant plasma formation occurs which sets in with an ionization avalanche. Note that the early time development including the ionization avalanche is obviously not influenced by the ionic correlations. Their development requires a significant amount of charged plasma particles which are only available after the ionization avalanche. The production of initial “seed charges” by black-body radiation does not depend on ionic correlations. The ionization avalanche appears somewhat ($`1\mu `$s) earlier than in the experiment since in reality the first electrons produced by the black-body radiation will leave the atom cloud until a sufficiently strong positive space charge builds up to trap the electrons. With the increasing number of free charges, the plasma forms quickly, followed by a partial back evolution into a Rydberg gas. Fig. 1 shows that this part of the gas dynamics may indeed be affected by ion-ion correlation, which leads to a more efficient recombination in the final stage of the plasma evolution. The reason for this behavior can be understood from eqs. (4) and (5). Together, $`W_c`$ and $`U_c`$ lead to an additional acceleration ($`2/9U_c`$ if Eq. (2) is used for $`u`$) in addition to the ideal thermal electron pressure $`k_\mathrm{B}T_e`$. Hence, the adiabatic electron cooling during the expansion is faster than without correlations. In turn recombination is accelerated, since the corresponding rate strongly increases with decreasing temperature. We note that corresponding effects of the same order of magnitude are also found for the experiments Kulin et al. (2000); Killian et al. (2001). The reason that they do not significantly affect the expansion dynamics of the plasma part lies in the fact that the correlations mainly influence the recombination into high-lying states at later stages of the evolution. Since these are very weakly bound, they hardly influence the kinetic energy of the system which determines the expansion velocity. Hence, good agreement is found between calculations with and without inclusion of correlation effects as long as quantities related to the plasma expansion are compared.

Figure 2 illustrates the dependence of the plasma evolution on the initial conditions. For fixed initial width of the Rydberg atom cloud (Fig. 2a), the plasma formation occurs earlier with increasing density, due to the fact that ionization is faster and more efficient at higher densities. The later stages of the evolution, on the other hand, depend only weakly on the density. This may be attributed to the fact that in the evolution equation (4), $`\rho `$ enters only in the form $`N/M`$, i.e. essentially as the degree of ionization. For fixed initial density $`\rho `$, the ionization avalanche occurs at the same time, in accordance with the argument made previously (Fig. 2b). However, for smaller width, the plasma expands faster, leading ultimately to lower electron temperatures and increased recombination. Finally, if the number of atoms $`N_a`$ is kept fixed the density changes when the size of the atom cloud is varied. Hence, Fig. 2c can be interpreted as the combined effect seen in Figs. 2a,b: The ionization avalanche occurs at different times for each curve (due to different atomic densities) and the relaxation leads to different final degrees of ionization (due to the different size of the atomic cloud).

Figure 3 shows the non-ionized part of the gas, more precisely the population of Rydberg levels in the initial and final stage of the expansion, as well as at the time of the maximum degree of ionization. Initially, all atoms are prepared with a principal quantum number $`n_0=70`$. At later times, the interplay between ionization, recombination and exciting and deexciting collisions leads to a decrease of the average excitation. This decrease in excitation is the main source of energy that triggers the ionization and the plasma expansion, energy absorbed from the black-body radiation field is negligible in comparison. Finally, recombination repopulates higher levels (which “freeze out” due to the rapid decrease of electron temperature), but the peak of the distribution is still found at relatively low excitations. It would be very interesting to measure the final distribution of Rydberg atoms, which should provide a stringent test for theoretical models.

In summary, we have described the formation of a cold plasma from a highly excited Rydberg gas, using a simple model based on kinetic equations for the plasma evolution. Microscopic molecular-dynamics simulations are limited to much shorter timescales due to the large numerical effort required. We have estimated the maximum effect of ionic correlations by including them into the model assuming instant equilibration of the ions. In this respect the calculations without equilibration (no correlation, ionic temperature zero) and with instant equilibration should bracket the actual ionization yield of the Rydberg gas. Given the difference in the ionization yields with and without ionic correlation in fig. 1, more elaborate calculations are worthwhile although at present it is unclear whether the plasma evolution including correlation can be followed over sufficently long times, e.g., with P3M-codes.

###### Acknowledgements.

It is a pleasure to thank P. Gould, E. Eyler, T. Killian and S. Kulin for helpful comments. Financial support by DFG through SPP 1116 is also gratefully acknowledged.


format: text

text:
# Current Short–Range Tests of the Gravitational Inverse Square Law

## I Introduction

There are several motivations for the experimental study of the gravitational force at short distances. Most generally, gravity is poorly understood below one centimeter. Up until the year 2000 it was not known if gravity obeyed the inverse square law below distance scales as large as several millimeters. Second (and closely related) is the opportunity to explore a large parameter space for new forces in nature. Even with the recent advances in short–distance measurements, experimental limits allow for new forces at least a million times stronger than gravity acting over distances resolvable to the unaided eye. Third, due to developments in theoretical particle physics, short–range gravity experiments have become laboratories to test for specific new physics beyond the Standard Model. Many recent models attempting to unify gravity and the other fundamental forces in the same theoretical framework predict modifications of gravity in the range of about a millimeter. Among these, predicted signatures of extra spatial dimensions have been most responsible for stimulating new experiments.

In these experiments, new phenomena, including the influence of extra dimensions, are expected to manifest themselves as departures from the Newtonian inverse square law. A common way to parameterize these effects is with the Yukawa interaction. The potential energy due to gravity and an additional Yukawa force between test masses $`m_1`$ and $`m_2`$ is given by:

$$V=\frac{Gm_1m_2}{r_{12}}[1+\alpha \mathrm{exp}(r_{12}/\lambda )],$$
(1)
where $`G`$ is the Newtonian gravitational constant, $`r_{12}`$ is the distance between the test masses, $`\alpha `$ is the strength of the new interaction relative to gravity, and $`\lambda `$ is the range. The Yukawa potential arises when the interaction between test masses $`m_1`$ and $`m_2`$ in the above expression is mediated by a massive field. In the case of a massless mediator with infinite range, such as the graviton, the usual $`1/r`$ potential is recovered. For cases in which the mediator couples to some quantity other than mass (for example, baryon number), $`\alpha `$ can acquire a dependence on the composition of the test masses. This effect is expected to be at most a small fraction of $`\alpha `$. It is addressed in experimental tests of the equivalence principle, which are conducted over ranges in which gravity is relatively well–measured.

This review concentrates on the search for new effects in the interaction range between a few microns and a few centimeters. As the experiments in this range are at or near gravitational sensitivity, it is here where phenomena associated with the extradimensional models of greatest interest are most likely to appear in the near future. There are also other specific predictions of new physics in this regime, some indirectly related to the physics of extra dimensions, which may soon be observable. Experimental programs motivated in part by the physics of extra dimensions have begun to explore other distance ranges as well. A significant effort is active in the range below 100 nm in connection with recent Casimir force measurements, as reviewed in Refs. fischbach  and bordag\_pr .

To date, no convincing signal of a departure from Newtonian gravity has been observed. In the absence of signals, experiments commonly report results in terms of limits on the strength parameter $`\alpha `$ for a given range $`\lambda `$. The current experimental limits are shown in Fig. 1, together with recent theoretical predictions.

## II Theoretical Motivations

Several predictions of new physics which can be tested in current short–range experiments arise in models addressing the unification of the fundamental forces. The leading candidate models are inspired by string or M-theory, which has to be formulated in more than three spatial dimensions. The extra dimensions are typically compactified at a scale on the order of the Planck length ($`L_P10^{35}`$ m or $`L_P^1=M_P10^{19}`$ GeV), the presumed gravitational unification scale, or concealed from direct detection by other mechanisms. Recently several models have been discovered which hinge on at least one of the extra dimensions remaining larger than $`L_P`$, with significant consequences for short–range experiments. String theories also contain scalar fields called moduli, which parameterize the size and shape of the extra dimensions and can mediate new short-range forces. Theoretical models are treated extensively in other parts of this volume; in this section we provide a brief description of some of the predictions.

In the recent model proposed by Arkani-Hamed, Dimopoulos, and Dvali (ADD), unification can occur at a scale $`L^1=M^{}<<M_P`$, with the relative weakness of the gravitational interaction arising as a consequence of $`n`$ “large” extra dimensions in which gravitons propagate but the Standard Model fields do not ADD . The size $`R`$ of the extra dimensions is related to the unification scale and the usual Planck scale by:

$$R^n=M_P^2/M^{(2+n)}$$
(2)
Setting $`M^{}=1`$ TeV, and effectively eliminating the hierarchy problem (the vast discrepancy or “desert” between the weak and Planck scales), the choice of $`n=1`$ implies $`R10^{13}`$ m, clearly ruled out by astrophysics. However, the choice of $`n=2`$ implies $`R1`$ mm, with the consequence that the gravitational force will behave according to a $`1/r^4`$ law below this scale. As the scale $`R`$ is approached from above, Yukawa corrections are predicted; the model illustrated in Fig. 1 predicts $`\alpha =4`$ floratos ; kehagias . In the absence of signals for these large compact dimensions, short–range gravity experiments can set limits on the unification scale $`M^{}`$, which is also probed directly in collider experiments.

A consequence of eliminating the hierarchy problem in favor of embedding the Standard Model on a 3-dimensional brane in a higher dimensional bulk is that the success of the “desert” in explaining the approximate symmetries responsible for a variety of phenomena (e.g., small neutrino masses and proton stability) no longer applies. In response to these issues, two authors of the ADD model consider a scenario in which the appropriate symmetry–breaking processes occur at a higher order on distant branes, with suppressed information about this breaking transmitted to the Standard Model brane via messenger fields arkani . The messenger fields should be easily detectable by a short–range macroscopic force search if they have mass of $`10^3`$ eV or less, as they are predicted to have couplings on the order of $`\alpha 10^6`$.

Several predictions arise in supersymmetric models with extra dimensions compactified on the order of the weak scale. Antoniadis, Dimopoulos and Dvali consider a model in which supersymmetry is broken at $`1`$ TeV by a compactification process, and show that the modulus associated with the dimension left at this scale acquires mass in the sub–eV range and a coupling of $`\alpha 1/3`$ antoniadis . The prediction for this “radius” modulus is shown in Fig. 1. Considering the problem of the stabilization of the extra dimensions in related models, Chacko and Perazzi predict a scalar field or “radion” corresponding to the size fluctuations of TeV–sized extra dimensions with a coupling of $`\alpha 30`$ and a range $`\lambda 40\mu `$chacko .

The other moduli predictions in Fig. 1 derive from an earlier model in which supersymmetry is broken near 100 TeV by a gauge–mediated process. This leads to moduli with masses in the sub–eV range and couplings as strong as $`\alpha 10^6`$ dimopoulos .

The dilaton is another scalar in string theory whose vacuum expectation value determines the string coupling constant. Dilaton mass can also arise from supersymmetry breaking but is generally sensitive to unknown physics. Predictions for the dilaton coupling vary widely in the literature; the lower and upper bounds in Fig. 1 are representative taylor ; ellis ; kaplan .

Fig. 1 also shows two predictions of new macroscopic effects not directly related to string theory. The axion, a light pseudoscalar motivated by the strong CP problem of QCD, is predicted to mediate forces between unpolarized test masses moody . The allowed range between 200 $`\mu `$m and 20 cm in Fig. 1 is constrained by laboratory and astrophysical observations rosenberg , and the coupling by measurements of the neutron electric dipole moment. The remaining prediction is motivated by the cosmological constant ($`\mathrm{\Lambda }`$) problem, the discrepancy between the observed flatness of the universe and the extreme curvature expected from the vacuum energy contributions of the Standard Model fields. Several models sundrum ; beane ; schmidhuber  have attempted to address this problem by introducing new interacting quanta with a range on the order of $`\mathrm{\Lambda }^{1/4}100\mu `$m.

## III Experimental Challenges

The principal limitation to testing gravity at short distances is the scaling of the signal force with the size of the apparatus. If all dimensions of a gravitational experiment are scaled by the same factor $`r`$, gravitational attraction behaves as $`r^4`$, so that signal forces become extremely weak at small distances. For example, expressing the gravitational force between two spherical test objects (radius $``$ separation $`r`$) in terms of their densities $`\rho `$ yields: $`FG\rho _1\rho _2r^4`$. For typical test mass densities ($`\rho 20\text{g/cm}^3`$), substituting $`r=10`$ cm gives $`F10^5`$ N. Measuring a force of this magnitude was the challenge met by Cavendish with his torsion balance experiment, published in 1798. In order to measure gravity at 100 $`\mu `$m, the experimenter is faced with scaling the size of the test masses to this range as well, so as not to be overwhelmed by the attraction of the additional mass at larger scales. Setting $`r=100\mu `$m yields $`F10^{17}`$ N, which explains in large part why it has taken two centuries to achieve measurements sensitive to gravitational effects below a millimeter.

At the same time, backgrounds due to known physics increase rapidly at short distances. Electrostatic effects due to surface potentials increase as $`r^2`$, and magnetic forces due to contaminants increase as $`r^4`$. At the lower end of the range considered ($`1\mu `$m), the Casimir force (of experimental interest itself) also becomes visible and quickly begins to dominate. The following sections are devoted to a description of current short–range experiments, with emphasis on the strategies used to overcome the dual challenge of extremely small signals and rapidly increasing backgrounds.

## IV Classical Gravitation Experiments

### The Irvine 2–5 cm Experiment and Cryogenic Torsion Pendulum

The highly linear and extremely sensitive torque response of thin filaments under tension has made the torsion balance the instrument of choice for measuring weak forces at laboratory scales. Between 3 mm and 3 cm, the best constraints on new forces derive from the 2–5 cm experiment in the group of R. Newman at the University of California at Irvine irvine . In this experiment, a cm–sized cylindrical test mass was suspended from one end of a torsion balance inside a long massive tube (a null geometry for $`1/r^2`$ forces). As the tube was translated periodically perpendicular to its axis, no deflection of the balance above that associated with a small Newtonian edge effect was observed. Torque sensitivities of about $`10^{14}`$ Nm were obtained, yielding the limit curve in Fig. 1.

A principal limitation to the sensitivity of this experiment was the effect of tilt. Motion of the tube along its guide track caused a slight tilt in the large block supporting the balance and a resulting phase–correlated motion of the balance suspension point. While the tilt was largely compensated with external instrumentation, the residual effect remained one of the main sources of systematic error. Magnetic effects presented another large source of systematic uncertainty, while torque measurement sensitivity was limited by a combination of instrument noise, seismic noise, and spurious gravitational forces not associated with the motion of the tube.

For the past several years, the Irvine group has been constructing a torsion pendulum which operates at cryogenic temperatures ($`T2`$ K), in collaboration with the group of P. Boynton at the University of Washington (Fig. 2). This approach has several advantages cryo , of which we mention only a few. For one, the pendulum may be operated in the frequency mode, which is highly insensitive to the effect of tilt. In this method, pioneered by the Boynton group Boynton , the resonant frequency of the pendulum is monitored for the shifts induced by the force gradient generated by an external source mass. At room temperature, the temperature dependence of the torsion fiber elastic properties limits the usefulness of this mode, but at cryogenic temperatures this dependence is small and temperature stability is easy to maintain. Operation at cryogenic temperature also allows for the use of superconducting shielding to reduce magnetic backgrounds.

Finally, cryogenic operation can reduce instrument noise considerably, inasmuch as this effect is dominated by thermal noise due to dissipation in the force–sensitive detector. In the thermal noise limit, the sensitivity of this experiment (and all others discussed below) to new forces is given by:

$$\alpha \frac{1}{\rho _1\rho _2}\sqrt{\frac{\omega T}{Q\tau }},$$
(3)
where $`\rho _1`$ and $`\rho _2`$ are the test mass densities, $`\omega `$ is the operational frequency of the experiment, $`T`$ is the temperature, $`Q`$ is the detector mechanical quality factor, and $`\tau `$ is the experimental integration time. In addition to the improvement in sensitivity expected from a reduction in $`T`$ from 300 K to 2 K, further improvement can be expected from a gain in $`Q`$ which is observed in many materials at low temperatures.

Construction of the cryogenic experiment is complete and a precision measurement of the Newtonian constant $`G`$ is in progress. A test of the inverse square law is planned which will have maximum sensitivity for forces with a range of about 15 cm. For the latter experiment, additional sensitivity is expected from a re-optimized test mass design based on a harmonic expansion of the mass and field multipole moments. In a method similar to that used in previous equivalence principle tests, the pendulum and source mass are configured to suppress Newtonian mass and field moments, respectively, so that classical gravitational effects enter into the measurement only to second order or higher moore . At the same time, non–Newtonian moments can be enhanced. If all noise sources can be reduced to near the level of the thermal noise, an improvement in sensitivity of about two orders of magnitude over previous results is expected, corresponding to the projected limit shown in Fig. 8.

### The Eöt-Wash Short-Range Experiment

Between $`100\mu `$m and 3 mm, the most sensitive limit on new forces has been attained by the short–range torsion pendulum experiment in the Eöt-Wash group at the University of Washington (Fig. 3washington\_prl . The force–sensitive part of the torsion pendulum consists of a 1 mm thick aluminum annulus with an array of ten equally spaced holes. The source mass consists of a stack of two copper disks, each a few millimeters thick, with similar arrays of holes. The source rotates approximately once every two hours, torquing the pendulum 10 times per revolution. This arrangement has several clever features. For one, the signal occurs at 10 times the revolution rate, allowing easy discrimination from vibrations associated with the source drive. Also, in the source stack, the angular position of the lower disk is offset relative to the upper disk so as to cancel the torque associated with the long–range Newtonian attraction between the test masses, while leaving possible sub–millimeter forces unaffected.

Not shown in Fig. 3, an essential component of the experiment is an electrostatic shield. This is made from a $`20\mu `$m thick beryllium–copper foil stretched on a frame between the source and pendulum.

The torque as a function of test mass separation was measured for two prototypes of this apparatus. Torque sensitivities of about $`10^{16}`$ Nm were attained, and no deviation from the expected Newtonian signal was observed for test mass separations down to about $`200\mu `$m. From these data, the preliminary results shown in Fig. 1 were obtained washington\_xxx . This is the shortest–ranged experiment to date to have attained gravitational sensitivity.

Improvements to this experiment currently under way include thinner pendulum and source disks, both constructed of copper or molybdenum for greater test mass density. The hole arrangements have been further optimized for Yukawa signals and cancellation of Newtonian gravity. Better vibration isolation and a thinner shield should significantly reduce the test mass separation. The projected limits expected with this technique are shown in Fig. 8.

## V High Frequency Techniques

With low operating frequencies and high mechanical quality factors ($`10^5`$), torsion balance experiments are very attractive from the point of view of thermal noise, even at room temperature. However, their sensitivities are usually limited by other backgrounds, including low–frequency vibrations which can also limit the minimum practical test mass separation. Recently, several experiments have been developed using high–frequency techniques, which show promise for operation at the thermal noise limit and for attaining smaller test mass separation.

### Colorado Torsional Oscillator Experiment

An approach for improving the experimental sensitivity at $`100\mu `$m and below has been pursued for the past several years by the authors and their collaborators at the University of Colorado colorado\_npb ; colorado\_mg9 . The basic design is shown in Fig. 4. Planar test mass geometry is used in order to concentrate as much of the available test mass density as possible at the range of interest. The flat plates are nominally a null geometry with respect to $`1/r^2`$ forces, which is important for suppressing Newtonian forces relative to new short–range effects. The test masses are fabricated from tungsten wafers, approximately 250 microns thick. The source mass consists of a 4 cm long cantilever, which is driven mechanically at a natural resonant frequency of the detector mass, nominally 1 kHz. The detector consists of two coplanar rectangles joined along their central axes by a short segment. In the resonant mode of interest, the detector rectangles counter–rotate in a torsional mode about the axis defined by the segment.

Electrostatic and acoustic backgrounds are suppressed with a stiff conducting shield, consisting of a gold–plated sapphire wafer, between the test masses. The shield currently limits the minimum test mass separation to approximately 100 microns.

Operation of the experiment at the detector resonant frequency, while increasing the experimental sensitivity, places strict requirements on vibration isolation. At the operational frequency of 1 kHz it is possible to construct a simple passive vibration isolation system, in which the test masses and shield are each suspended from a series of brass disks connected by fine wires under tension. This system provides about 200 dB of attenuation at the operational frequency chan . Modular design, visibility, and easy control of test mass positioning and electrical connections allow for rapid characterization of observed signals in situ.

Electromagnetic, acoustic and vibrational backgrounds have been sufficiently suppressed such that the experiment is currently limited by detector thermal noise. To reduce the thermal background, the detector is annealed at 1300 C for several hours before installation in the experiment, which significantly increases its mechanical quality factor.

In a recent series of runs, no signal above thermal noise was observed after 22 hours of integration time. This corresponded to a force sensitivity of about $`10^{15}`$ N, from which the limit derived in Fig. 1 is obtained colorado\_xxx . Development of a stretched–membrane shield and flatter test masses with higher Q is in progress. In the absence of backgrounds, the limits attainable with the current sensitivity of the experiment (at 300 K) at half the present minimum test mass separation are shown in Fig. 8. A cryogenic version of this experiment, with a 4 K operating temperature and an expected detector $`Q`$ enhancement of an order of magnitude, could improve the sensitivity by an additional factor of 40.

### Stanford Micro-cantilever Experiment

Another high–frequency measurement is in progress in the group of A. Kapitulnik at Stanford. This experiment uses techniques of scanning force microscopy and microfabrication to probe distance scales on the order of a few tens of microns (Fig. 5). The force–sensitive detector mass consists of a $`50\mu `$m gold cube mounted on the end of a silicon nitride cantilever $`250\mu `$m long and less than $`1\mu `$m thick.

The source mass consists of a parallel array of alternating strips of gold and silicon, each $`100\mu `$m wide. It is scanned in the horizontal plane below the detector in the direction perpendicular to the strips, with the resulting periodic variations in density providing the driving mass–coupled force. This design offers an advantage similar to that of the Eöt-Wash torsion pendulum, in that the frequency of the source drive mechanism can be set much lower than that of the expected signal. In the Stanford experiment, the mass modulation frequency is matched to the detector cantilever resonance (about 300 Hz), resulting in less of a burden on vibration isolation than for gap–modulated experiments.

An electrostatic shield, fabricated from a silicon nitride membrane $`3\mu `$m thick with a 100 nm gold plating, is placed between the test masses. It is attached to the base of the cantilever mount, leaving a gap of $`15\mu `$m below the detector mass.

The limiting background to this measurement is expected to be thermal noise of the force–sensitive cantilever. To reduce this background, the experiment is inserted into a liquid helium dewar for operation at $`T10`$ K, where the expected force sensitivity is about $`10^{17}`$ N for integration times of a few hours.

After 20 minutes of integration time during a run in 2002, data from this experiment showed evidence of a signal approximately three times greater than the expected thermal noise stanford . The signal appeared most likely to be an electrostatic background, due to modulation of the shield caused by surface variations on the source, and a 0.3 V potential between the shield and the detector mass. The phase and magnitude of the observed signal was not consistent with a mass–coupled force generated by the source mass, and the data were used to obtain the preliminary limits shown in Fig. 1 representing an upper limit on such forces. While a projected limit was not available as of this writing, improvements to the experiment continue, including a new source mass with additional conducting and insulating layers above the strips to eliminate surface variations.

### Dusseldorf Torsional Oscillator Experiment

An experiment in progress in the group of S. Schiller at the University of Dusseldorf introduces several new features. An early conception is illustrated in Fig. 6. Planar geometry is used, with both test masses mounted vertically. The detector mass, fabricated from a high–purity silicon wafer, is a torsional oscillator similar in design to the Colorado detector. Excitation of the source mass at the detector resonance (about 5 kHz) is accomplished using a rotating disk with nine tungsten alloy teeth embedded in a lower–density substrate.

A thin metal plate is placed between the test masses as an electrostatic shield. All elements are mounted atop multi–stage vibration isolation stacks, consisting of viton rings sandwiched between 1 kg–mass brass disks.

The sensitivity of this experiment is also expected to be limited by thermal noise. In the thermal–noise limit, the reduction in sensitivity associated with the use of a low–density silicon detector is offset somewhat by the exceptionally high mechanical quality factors attainable, which are about an order of magnitude larger than those observed so far in metal oscillators at room temperature. Plans for operation of the experiment at 4 K are under way, where improvement in $`Q`$ of at least two additional orders of magnitude can be expected for silicon detectors. The combined effect of the increase in $`Q`$ with the reduction in temperature is reflected in the very sensitive projected limits for this experiment, shown in Fig. 8.

A room–temperature prototype of the experiment has been assembled and is under test. In early runs, a signal approximately 10 times larger than the detector thermal noise was observed. As the signal was present even for large separations of the test masses, a vibrational background was suspected. This background has been reduced with several improvements in the design, including a more uniform source disk and remote positioning of the drive vibration isolation stack. A cryostat for the 4 K version of the experiment is also under construction.

### Padova Casimir Force Experiment

Below test mass separations of a few microns, the principal background to macroscopic force searches is expected to be the Casimir effect, a force which arises due to zero–point fluctuations of the electromagnetic field. At the moment, the best limits on new forces in this range are derived from direct measurements of this effect, in general by ascertaining the maximum possible size of a Yukawa–type force for the Casimir plates still consistent with the experimental residuals. Below $`\lambda =3\mu `$m, the best published limits are obtained from the measurement by Lamoreaux in 1997 bordag\_pr ; colorado\_npb  which used a torsion balance to measure the Casimir force between a flat gold–plated disk and a thin spherical shell Lamoreaux .

An experiment at Legnaro laboratories by a joint team of the INFN, Padova and Pavia, and the University of Padova, uses a frequency–shift technique to measure the Casimir force between parallel conducting plates (Fig. 7padova\_prl . The force–sensitive detector plate, a $`47\mu `$m thick, 2 cm long silicon reed with a 50 nm chromium plating, is clamped at its base and driven electrically into oscillation at its lowest natural frequency of about 138 Hz. The source mass, a 5 mm thick, chromium–plated silicon block, is mounted on a nearby translation stage and brought into very close proximity to the detector. The detector resonance is monitored optically for the frequency shifts induced by static external force gradients.

The apparatus is mounted on passive vibration isolation stacks, and in turn into a large scanning electron microscope (not shown in Fig. 7) in order to monitor the cleanliness and parallelism of the test mass surfaces at the micron–sized separations needed for this measurement. The electrostatic background is due primarily to contact potentials in the circuit connecting the test masses. It is eliminated to first order by nulling the static force between the test masses at small distances, and the remaining background is subtracted off after careful fits to the measured force in the presence of applied potentials at larger gaps.

This experiment was able to obtain a measurement of the Casimir force to a precision of 15%, the first unambiguous measurement of this effect for the parallel–plate geometry. An earlier version of this experiment was among the first high–frequency measurements used to derive limits on new effects padova\_prd . While those limits were not competitive with the most stringent results at that time, a re–optimized version of the present apparatus is expected to improve the limits by at least an order of magnitude in the range near $`1\mu `$m, as shown in Fig. 8 padova\_cqg .

## VI Conclusion

Motivated in large part by suggestions of extra dimensions, three of the experiments described above have explored approximately 7 square decades of the parameter space for new forces in the past five years historical , and have set significant limits on new physics. The maximum possible size of any extra dimension in the ADD scenario has been limited to $`200\mu `$m (or $`150\mu `$m for $`n=2`$ equal–sized extra dimensions, corresponding to a lower limit on the unification scale of $`M^{}4`$ TeV) washington\_xxx , a new upper limit on the dilaton mass has been set at $`8.6\times 10^3`$ eV (assuming a limit on the coupling of $`\alpha 2000`$), and much of the remaining parameter space for the moduli force scenario of Ref. dimopoulos  has been excluded.

While the torsion pendulum remains the instrument of choice for scales above 100 microns, high–frequency experiments show great potential for improving the limits at shorter ranges. Fig. 8 shows the projected limits for the next few years, considering only some of the experiments currently in progress. Experimental coverage of at least 7 additional square decades of parameter space, with investigation of the region of cosmological interest at gravitational strength, seems well within reach in the near future.

## Acknowledgements

We wish to thank R. Newman, P. Boynton, and M. Moore for description of the cryogenic torsion pendulum experiment, L. Haiberger and S. Schiller for materials describing the Dusseldorf experiment, and E. Adelberger, J. Chiaverini, C. D. Hoyle, A. Kapitulnik, R. Onofrio, and G. Ruoso for useful materials and comments.


text:
# The Large 𝑁 Glueball Mass Spectrum in 2+1 Dimensions

## I Introduction

In this paper we study the glueball mass spectrum for large $`N`$ pure gauge theory (no quarks) in 2+1 dimensions using Hamiltonian LGT. Our primary motivation lies in recent developments in string theory where a direct connection between pure gauge theory in the large $`N`$ limit and certain string theories has been conjectured. An independent test of this conjecture is provided by large $`N`$ calculations of the glueball mass spectrum in LGT. Impressive Monte Carlo calculations by Teper Teper (1999) and Lucini and Teper Lucini and Teper (2002, 2001) have progressed to the stage where accurate values of glueball masses have been calculated up to $`N=6`$ in 2+1 dimensions. This has allowed an extrapolation to the $`N\mathrm{}`$ limit. With the analytic techniques of the Hamiltonian approach Carlsson and McKellar (2003) it is possible to extend these calculations to much larger $`N`$, allowing a more reliable $`N\mathrm{}`$ extrapolation.

The outline of this paper is as follows. We start in Section II with a historical background, where the original motivations for large $`N`$ physics are discussed and the recent developments in string theory introduced. In Section III we use the analytic variational approach discussed in Ref. Carlsson and McKellar, 2003 to calculate glueball masses at finite values of $`N`$ up to $`N=25`$ in some cases. These results are then extrapolated to explore the large $`N`$ limit of the 2+1 dimensional glueball mass spectrum. The $`N\mathrm{}`$ limits obtained here are compared to the results of the analogous Monte Carlo studies in the Lagrangian approach.

## II Background

The nonperturbative region of quantum chromodynamics (QCD) has been examined predominantly with numerical techniques. While the numerical Monte Carlo simulations of Lagrangian LGT have made significant progress, there have been few developments in analytic techniques. One analytic technique which has received considerable attention in recent years is the large $`N`$ limit.

In 1974 ’t Hooft proposed a study of SU($`N`$) gauge theories in the large $`N`$ limit rather than the physically interesting $`N=3`$ case ’t Hooft (1974). It was hoped that the SU($`N`$) theory could be solved analytically and would be, in some sense, close to SU(3). Indeed it was shown that the large $`N`$ theory simplifies drastically ’t Hooft (1974) but still captures at least some of the complexity of the SU(3) theory Witten (1979); Eguchi and Kawai (1982). Motivated by the fact that the QCD coupling constant $`g^2`$ is a poor expansion parameter because its value depends on the energy scale of the process under consideration, t’ Hooft considered an expansion in another dimensionless parameter of SU($`N`$) QCD, $`1/N`$. Based on an ingenious organisation of Feynman diagrams, ’t Hooft was able to show that in the large $`N`$ limit, keeping the ’t Hooft coupling $`g^2N`$ fixed, only planar diagrams remain. While the remaining planar theory can be solved exactly in two dimensions, the three and four dimensional theories have not yet been solved.

Recent developments in string theories have attempted to address this problem. In 1997 Maldacena conjectured a correspondence, in the duality sense, between superconformal field theories and string theory propagating in a non-trivial geometry Maldacena (1998). The idea behind duality is that a single theory may have two (or more) descriptions with the property that when one is strongly coupled the other is weakly coupled. A technique for breaking the conformal invariance and supersymmetry constraints of Maldacena’s original proposal was later developed by Witten Witten (1998). This led to the hope that nonperturbative pure SU($`N`$) gauge theories could be described analytically in 3+1 dimensions by their string theory dual. To develop this hope into a concrete solution would require perturbative expansions within the dual string theory. Such developments would seem to be some way off.

In order to test the viability of using string theories to probe nonperturbative QCD independent tests are needed. The glueball mass spectrum of QCD provides a perfect laboratory. Recent Monte Carlo calculations have provided stable estimates of glueball masses up to $`N=6`$ Teper (1999); Lucini and Teper (2001) in 2+1 dimensions and up to $`N=5`$ Lucini and Teper (2001) in 3+1 dimensions. From these estimates an $`N\mathrm{}`$ limit can be extrapolated in each case. Convincing evidence of $`𝒪(1/N^2)`$ finite $`N`$ corrections have been demonstrated in 2+1 dimensions, verifying a specific prediction of ’t Hooft’s $`1/N`$ expansion in the quarkless case. From the string theory side, present calculations require a strong coupling limit ($`N,g^2N\mathrm{}`$) to be taken. In this limit the string theory reduces to classical supergravity and the calculation of glueball masses, $`m(J^{PC})`$, is straightforward. In this limit the glueball is represented by a dilaton field whose mass can be extracted by solving the dilaton wave equation Csaki et al. (1999a); de Mello Koch et al. (1998); Zyskin (1998). Despite the extreme approximations required surprising quantitative agreement with the weak coupling results from LGT were obtained. Later it was realised that, in 2+1 dimensions, lower mass states can be constructed from the graviton field Brower et al. (2000a). With this approach the quantitative agreement with LGT is lost but qualitative agreement in the form of the prediction

$`m(0^{++})<m(2^{++})<m(1^+)`$ (1)
still holds. The qualitative agreement between string theories and LGT persists in 3+1 dimensions Brower et al. (2000b) with

$`m(0^{++})<m(2^{++})<m(0^+),`$ (2)
for both theories. It appears that the identification of glueball states remains a problem. A serious difficulty is the removal of unwanted states which have masses of the same magnitude as the glueballs. Recipes for removing these unwanted states have been proposed Russo (1999); Csaki et al. (1999b) but new parameters need to be introduced.

The exploration of large $`N`$ gauge theories outside of string theory has not been restricted to traditional Monte Carlo methods. Recent studies using the light-front Hamiltonian of transverse LGT have shown agreement with the Monte Carlo calculations of Teper Dalley and van de Sande (2001). In this approach explicit calculations at $`N\mathrm{}`$ are possible without the need for extrapolation. Continuum calculations, without the use of string theory, have also commenced. An impressive series of papers Karabali and Nair (1996); Karabali et al. (1998a, b) has led to specific predictions for the string tension for all $`N`$ in 2+1 dimensions Karabali et al. (1998b). These predictions lie within $`3\%`$ of lattice calculations up to $`N=6`$ beyond which LGT results are not available. A recent review is available in Ref. Nair, 2002.

## III Extrapolation

### III.1 Introduction

It this section we calculate 2+1 dimensional SU($`N`$) glueball masses for various $`N`$ using the analytic variational technique of Ref. Carlsson and McKellar, 2003. Our intention is to explore the large $`N`$ limit of glueball mass spectrum by extrapolating the results obtained to their $`N\mathrm{}`$ limit. Employing analytic methods we can extend the calculation of the lowest energy massgaps to SU(25) on a desktop computer. This should be compared to current results in the Lagrangian approach which are limited, at present, to SU(6). With large values of $`N`$ available a more reliable extrapolation to the $`N\mathrm{}`$ limit is possible.

To simulate the ground, or (perturbed) vacuum, state with energy $`E_0`$, we use the one plaquette trial state,

$$\begin{array}{c}\text{}\end{array}$$
(3)
Here, $`|0`$ is the strong coupling vacuum defined by $`_i^\alpha (𝒙)|0=0`$ for all $`i`$, $`𝒙`$ and $`\alpha =1,2,\mathrm{},N^21`$. $`_i^\alpha (𝒙)`$ is the lattice chromoelectric field on the directed link running from the lattice site labelled by $`𝒙`$ to that labelled by $`𝒙+a\widehat{𝒊}`$. Here $`\widehat{𝒊}`$ is a unit vector in the $`i`$ direction. The directed square denotes the traced ordered product of link operators, $`U_i(𝒙)`$, around an elementary square, or plaquette, of the lattice,

$$\begin{array}{c}\text{}\end{array}$$
(4)
where $`a`$ is the lattice spacing.

As discussed in Ref. Carlsson and McKellar, 2003, the precise form of $`c(\beta )`$, where $`\beta =N/g^2`$, is unimportant in a scaling region. As it turns out, to calculate variational wave functions, by minimising the vacuum energy density, for large dimension gauge groups is cumbersome. Numerical precision becomes a factor in the minimisation of the energy density. This problem is magnified in the calculation of tadpole improved results using the iterative procedure employed in Ref. Carlsson and McKellar, 2003. We thus abandon the use of variational wave functions here. Instead we make use of the one plaquette wave function of Eq. (3) and define a simple dependence of $`c(\beta )`$. Here we choose $`c(\beta )=\beta `$.

Many Hamiltonians are available for LGT calculations Carlsson and McKellar (2001). In Ref. Carlsson and McKellar, 2003 we used Kogut-Susskind Kogut and Susskind (1975), improved and tadpole improved Hamiltonians to calculate glueball masses for $`N5`$. Here we employ the simplest of these, the Kogut-Susskind Hamiltonian, which is defined for pure SU($`N`$) gauge theory with coupling $`g^2`$ on a lattice with spacing $`a`$ by

$``$ $`=`$ $`{\displaystyle \frac{g^2}{2a}}{\displaystyle \underset{𝒙,i}{}}_i^\alpha (𝒙)^2+{\displaystyle \frac{2N}{ag^2}}{\displaystyle \underset{𝒙,i<j}{}}P_{ij}(𝒙),`$ (5)
where the plaquette operator is given by

$$\begin{array}{c}\text{}\end{array}$$
(6)
To evaluate the matrix elements appearing in the glueball mass calculation we use the SU($`N`$) generating functions derived in Ref. Carlsson and McKellar, 2003,

$`G_{\mathrm{SU}(N)}(c,d)`$ $`=`$ $`{\displaystyle _{\mathrm{SU}(N)}}𝑑Ue^{c\mathrm{Tr}U+d\mathrm{Tr}U^{}}`$ (7)
$`=`$ $`{\displaystyle \underset{l=\mathrm{}}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{d}{c}}\right)^{lN/2}det\left[I_{l+ji}\left(2\sqrt{cd}\right)\right]_{ii,jN}`$
and

$`H_m(c,d)`$ $`=`$ $`{\displaystyle _{\mathrm{SU}(N)}}𝑑U\mathrm{exp}\left[c(\mathrm{Tr}U+\mathrm{Tr}U^{})+d\mathrm{Tr}(U^m)\right]`$ (8)
$`=`$ $`{\displaystyle \underset{l=\mathrm{}}{\overset{\mathrm{}}{}}}det\left[{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{d^k}{k!}}I_{l+ji+mk}(2c)\right]_{1i,jN}.`$
Here the quantities inside the determinant are to be interpreted as the $`(i,j)`$-th entry of an $`N\times N`$ matrix.

For the minimisation process, we use the same basis of rectangular states as in Ref. Carlsson and McKellar, 2003. With this basis the states $`J^{++}`$ and $`J^{}`$, with $`J=0`$ or 2, are accessible. The $`J=2`$ states require the use of a minimisation basis that contains no square states so that the excited states are invariant under rotations by $`\pi `$ but not $`\pi /2`$. The dependence of the massgaps on the coupling does not change significantly as $`N`$ is increased. We do however observe the appearance of what could possibly be additional low $`\beta `$ scaling regions as $`N`$ is increased.

### III.2 Convergence with $`𝒍_{\mathrm{𝐦𝐚𝐱}}`$

In Ref. Carlsson and McKellar, 2003 the dependence of the variational parameter on the truncation of the $`l`$-sum in Eq. (7) was considered. The calculation of the variational parameter depends only on the plaquette expectation value. The dependence of the massgaps on the truncation should also be explored as their calculation incorporates additional matrix elements. We find that fast convergence is achieved as $`l_{\mathrm{max}}`$ is increased for the $`0^{++}`$ massgaps. A typical example is shown in Fig. 1(a), where the massgaps are almost indistinguishable up to $`\beta =80`$ on the scale of the plot. For $`0^{}`$ the convergence is not as fast. A typical example is shown in Fig. 1(b). We find that in order to obtain convergence of the antisymmetric massgaps up to $`1/g^2350`$ for SU(25) we need a truncation of $`l_{\mathrm{max}}=30`$.

### III.3 The small $`𝜷`$ minima

The massgaps for $`N2`$, in units of $`Ng^2/a`$ do not differ significantly as functions of $`1/g^2`$. The lowest lying $`0^{++}`$ glueball masses have the characteristic form shown in Fig. 1(a) for each $`N`$. We observe a minimum at small $`\beta `$ and a scaling plateau at large $`\beta `$.

In this section we examine the minima occurring in the lowest $`0^{++}`$ and $`0^{}`$ glueball masses and consider the possibility that they may correspond to a scaling region. It appears that for SU(3) this region has been interpreted as a scaling region by some authors Fang et al. (2002); Chen et al. (1995); Leonard (2001). From the scaling arguments of Ref. Carlsson and McKellar, 2003, by changing the functional dependence of the parameter appearing in the vacuum state the scaling region can be modified. In this way it is possible to choose the variational parameter’s dependence on $`\beta `$ such that the region in which the mass gap takes its minimum is stretched over a large $`\beta `$ interval. In the same way, the scaling plateaux occurring for large $`\beta `$ may be compressed. With this in mind, it is unclear which scaling region we should take to be the correct one. For this reason we analyse both potential scaling regions. We start, in this section, with the small $`\beta `$ minima and consider the large $`\beta `$ plateaux in the next section.

Fig. 2 shows the lowest lying SU($`N`$) $`0^{++}`$ glueball masses in 2+1 dimensions in units of $`g^2N/a`$ as functions of $`1/g^2`$ for various $`N[3,25]`$. It is apparent that the minima depend only weakly on $`N`$. We see that the minima appear to approach a finite limit from below as $`N\mathrm{}`$. On the scale of Fig. 2 the minima corresponding to $`N=15`$ and $`N=25`$ are barely distinguishable.

Let us denote the minima in Fig. 2 by $`a\mathrm{\Delta }M^c/(Ng^2)`$ and consider them as a functions of $`N`$. Fitting $`a\mathrm{\Delta }M^c/(Ng^2)`$ to the model

$`{\displaystyle \frac{a\mathrm{\Delta }M^c}{Ng^2}}=\gamma _1+{\displaystyle \frac{\gamma _2}{N^{2\gamma _3}}}`$ (9)
for $`N8`$ gives best fit parameters

$`\gamma _1`$ $`=`$ $`0.83262\pm 0.00022`$
$`\gamma _2`$ $`=`$ $`0.97\pm 0.13`$
$`\gamma _3`$ $`=`$ $`0.990\pm 0.035.`$ (10)
We fit on the data $`N8`$ to minimise pollution from next to leading order corrections. If we assert that the power of the leading order correction for finite $`N`$ must be an integer power of $`1/N`$, then from Eq. (10) that order must be 2 in agreement with the predictions of large $`N`$ QCD. Fitting to a model with $`1/N^2`$ corrections for $`N8`$ gives

$`\mu _1(N^2)=0.83256\pm 0.00007{\displaystyle \frac{0.9753\pm 0.0072}{N^2}}.`$ (11)
We can go further and attempt to find the next to leading order corrections by fitting to the model

$`{\displaystyle \frac{a\mathrm{\Delta }M^c}{Ng^2}}=\gamma _1+{\displaystyle \frac{\gamma _2}{N^2}}+{\displaystyle \frac{\gamma _3}{N^{2\gamma _4}}}.`$ (12)
The best fit parameters when fit on the whole data set are

$`\gamma _1`$ $`=`$ $`0.83287\pm 0.00011`$ (13)
$`\gamma _2`$ $`=`$ $`1.116\pm 0.028`$ (14)
$`\gamma _3`$ $`=`$ $`1.715\pm 0.085`$ (15)
$`\gamma _4`$ $`=`$ $`1.627\pm 0.050`$ (16)
which is in disagreement with integer power next to leading order corrections. A good fit however is achieved to the model with $`1/N^4`$ next to leading order corrections:

$`\mu _2(N^2)`$ $`=`$ $`0.83233\pm 0.00018{\displaystyle \frac{0.986\pm 0.010}{N^2}}+{\displaystyle \frac{2.772\pm 0.090}{N^4}}.`$ (17)
The precise locations of the minima, $`a\mathrm{\Delta }M^c/(Ng^2)`$, are plotted as a funtion of $`1/N^2`$ in Fig. 4 along with the fitted models of Eqs. (11) and (17).

From Eq. (11), our best fit result for the $`N\mathrm{}`$ limit of $`a\mathrm{\Delta }M^c/(Ng^2)`$ is $`0.83256\pm 0.00007`$. The error here is purely statistical. We should, as always, expect a significant systematic error attributable to our choice of ground state and minimisation basis. This result may be compared to the Monte Carlo calculation of Lucini and Teper who obtain the $`N\mathrm{}`$ limit of the lowest $`0^{++}`$ glueball mass as $`0.8116\pm 0.0036`$ Lucini and Teper (2002). This is comparable to the result presented here. The result of Lucini and Teper is based on a linear extrapolation to the $`N\mathrm{}`$ limit of $`2N6`$ data. Their result differs significantly from the one presented here at the leading order finite $`N`$ corrections. The correction term obtained by Lucini and Teper, $`(0.090\pm 0.028)/N^2`$, has the same sign but is significantly smaller than the one presented here. When our data is fit on the range $`3N6`$ (we don’t obtain a small $`\beta `$ minimum for SU(2)) our slope is halved but is still significantly larger than that of Lucini and Teper. It should be pointed out that Lucini and Teper’s calculation was performed in the Lagrangian approach in which the coupling is the so called Euclidean coupling $`g_E`$. The Euclidean coupling and the Hamiltonian coupling, $`g^2`$, are equal up to order $`g_E^2`$ corrections. The precise relation between the couplings is only known for small $`g_E^2`$ Hamer (1996), a case which does not apply for the small $`\beta `$ minima. It appears that our simple calculation induces a level crossing whereby the lowest mass glueball state switches to a higher mass state beyond the small $`\beta `$ minimum. Although we do not have an explanation for this, presumably, by including additional states in the minimisation basis or implementing a more complicated vacuum state, the level crossing would no longer appear and the small $`\beta `$ minima would extend into large $`\beta `$ scaling regions. This should be checked in further studies. A glueball mass extracted from a large $`\beta `$ scaling region can be confidently compared to a corresponding Lagrangian calculation, for in that region of couplings the ratio of $`g_E^2`$ to $`g^2`$ is unity.

We obtain small $`\beta `$ minima for the $`0^{}`$ massgaps but only for $`N5`$. Plots of these massgaps with $`N5`$ are shown in Fig. 3. The behaviour is significantly different to that observed for the symmetric massgap. By $`N=25`$ the minima do not appear close to convergence. Indeed if convergence is occurring at all, it is significantly slower than was observed for the $`0^{++}`$ state. To explore this further, in Fig. 5 we plot the minima of Fig. 3, which we again denote by $`a\mathrm{\Delta }M^c/(Ng^2)`$, as a function of $`1/\sqrt{N}`$. We observe linear behaviour in the large $`N`$ limit. Fitting the $`N9`$ data to the model

$`{\displaystyle \frac{a\mathrm{\Delta }M^c}{Ng^2}}=\gamma _1+{\displaystyle \frac{\gamma _2}{N^{\gamma _3/2}}},`$ (18)
gives best fit parameters

$`\gamma _1`$ $`=`$ $`0.41\pm 0.01`$
$`\gamma _2`$ $`=`$ $`1.25\pm 0.02`$
$`\gamma _3`$ $`=`$ $`0.98\pm 0.04,`$ (19)
which is consistent with $`\gamma _3=1`$. Fitting the $`N9`$ data to a model with leading order finite $`N`$ corrections starting at $`1/\sqrt{N}`$ gives

$`\nu _1(N)=0.41390\pm 0.00007+{\displaystyle \frac{1.255\pm 0.003}{\sqrt{N}}}.`$ (20)
As was done for the symmetric case we can attempt to go further and find the form of the next to leading order finite $`N`$ corrections. Fitting the $`N5`$ data to the model

$`\nu _2(N)=\gamma _1+{\displaystyle \frac{\gamma _2}{\sqrt{N}}}+{\displaystyle \frac{\gamma _3}{N^{\gamma _4/2}}}.`$ (21)
gives best fit parameters

$`\gamma _1`$ $`=`$ $`0.410\pm 0.001`$ (22)
$`\gamma _2`$ $`=`$ $`1.2718\pm 0.0045`$ (23)
$`\gamma _3`$ $`=`$ $`76\pm 23`$ (24)
$`\gamma _4`$ $`=`$ $`9.5\pm 0.4`$ (25)
indicating vanishingly small next to leading order corrections for $`N5`$. The models $`\nu _1`$ and $`\nu _2`$ are plotted against the $`N5`$ data in Fig. 5.

The small $`\beta `$ minima appear to converge to a non-zero $`N\mathrm{}`$ limit. The convergence is significantly slower than for the small $`\beta `$ minima in the symmetric case. It is also clear that the leading order finite $`N`$ corrections are not the expected $`𝒪(1/N^2)`$ but rather $`𝒪(1/\sqrt{N})`$. The $`N\mathrm{}`$ limit of $`0.41\pm 0.01`$ is approximately half the corresponding $`0^{++}`$ result. Such a state does not appear in the calculation of Teper.

It is possible that the small $`\beta `$ minima are spurious. Dimensional analysis gives the expected scaling form but does not allow us to decide which scaling region is preferable, or if one is a lattice artifact. However, based on lattice calculations to date, we expect the $`0^{++}`$ glueball to be lighter than the $`0^{}`$ in the continuum limit. Spurious scaling regions have been observed in other lattice calculations Wichmann (2001). It is clear that additional analysis is needed here. Again, an important step will be to include additional non-rectangular states in the minimisation basis.

### III.4 The large $`𝜷`$ plateaux

Having considered the small $`\beta `$ minima as possible scaling regions we now move on to the large $`\beta `$ plateaux. This scaling region appears for all $`N`$ and for all states considered and so its interpretation as a genuine scaling region is less dubious. Furthermore, the glueball mass results extracted from these scaling regions may be confidently compared to the corresponding Lagrangian calculations since the ratio of the Hamiltonian to Euclidian coupling is unity up to small $`𝒪(g_E^2)`$ corrections.

We start with the $`0^{++}`$ state for which the best scaling behaviour is obtained. We calculate the five lowest lying massgaps corresponding to the five lowest glueball masses accessible with our choice of ground state and minimisation basis. We consider values of $`N`$ in the range $`3N25`$. For each $`N`$ considered we find a large $`\beta `$ scaling plateau for each of the five lowest mass states. The lowest lying massgap is shown in Fig. 6 for a range of $`N`$. We observe that in units of $`Ng^2/a`$ the massgaps do not depend strongly on $`N`$ and that in the scaling region they appear to approach a finite limit. Similar observations can be made for the four higher mass states obtained although the scalaing behaviour is less precise. We show the second lowest glueball mass in Fig. 7 as an example. A continuum limit for each massgap is obtained in the scaling region by fitting to a constant. For each fit we use a region of at least 10 data points which minimises the standard error. The continuum limit results for the five lowest lying massgaps, denoted by $`a\mathrm{\Delta }M^c/(Ng^2)`$, are shown as functions of $`1/N^2`$ in Figs. 6 and 9. Also shown in the plots are fits to models with leading order large $`N`$ corrections starting $`1/N^2`$,

$`\kappa _i^{++}=p_i+{\displaystyle \frac{q_i}{N^2}}+{\displaystyle \frac{r_i}{N^4}}.`$ (26)
Here the subscript $`i`$ labels the $`i`$-th lowest glueball mass. The values of the best fit parameters are given in Table 1.

These results should be compared with those of Lucini and Teper who obtain large $`N`$ glueball masses in units of $`Ng^2/a`$ in the $`0^{++}`$ sector with linear fits given by

$`0^{++}:`$ $`0.8116(36){\displaystyle \frac{0.090(28)}{N^2}}`$
$`0^{++}:`$ $`1.227(9){\displaystyle \frac{0.343(82)}{N^2}}`$
$`0^{++}:`$ $`1.65(4){\displaystyle \frac{2.2(7)}{N^2}}.`$ (27)
The fit for the $`0^{++}`$ state is obtained using $`4N6`$ data. The remaining fits are obtained using $`2N6`$ data. We find that the lowest glueball mass extracted from our large $`\beta `$ plateaux is consistent, in the $`N\mathrm{}`$ limit, with the state which Lucini and Teper label $`0^{++}`$. The slopes of the fits are of the same sign but differ significantly in magnitude. The $`0^{++}`$ state of Lucini and Teper does not appear in our data in the form of a large $`\beta `$ scaling plateau. There is however a hint of an approach to scaling in the vicinity of their result in our second lowest massgap data. This effect is only visible in our data for $`N13`$ and occurs for small $`\beta `$ as shown in Fig. 10. The values of $`\beta `$ for which this effect appears are quite close to those for which the small beta minima are observed in the lowest mass eigenstate. Similar effects are observed in the $`0^{}`$ data but the effect is much less convincing with the data currently available.

Let us consider this small $`\beta `$ region as a possible scaling region. We fit a constant to the available $`N15`$ data on a range of at least 4 data points chosen so that the standard error is minimised. These scaling values are plotted as a function of $`1/N^2`$ in Fig. 11. Also shown is the best fit linear model

$`\kappa _2^{++}=1.7605\pm 0.0032{\displaystyle \frac{5.83\pm 0.97}{N^2}}.`$ (28)
This produces an $`N\mathrm{}`$ limit which is close to but inconsistent with the result obtained by Lucini and Teper for their $`0^{++}`$ state.

We now move on to the $`0^{}`$ states. For these states the scaling behaviour is not as precise as that obtained for the $`0^{++}`$ states. Again we calculate the five lowest lying massgaps corresponding to the five lowest mass glueballs accessible with our choice of minimisation basis and ground state. For each $`N`$ considered up to 25 we find a large $`\beta `$ scaling plateau for each of the five states. The lowest lying massgap is shown in Fig. 12 for a range of $`N`$. The second lowest energy massgap is shown in Fig. 13. Considerably less data has been obtained for the $`0^{}`$ states due to the large $`l_{\mathrm{max}}`$ required for convergence as discussed in Section III.2. Despite this, the results qualitatively replicate those of the $`0^{++}`$ states. We observe that in units of $`Ng^2/a`$ the massgaps do not depend strongly on $`N`$ and that in the scaling region they appear to approach a finite limit. Similar observations can be made for the four higher mass states although the scaling behaviour is less precise. The continuum limit values extracted are plotted as functions of $`1/N^2`$ in Figs. 14 and 15. The dashed lines show fits of the continuum limit values to the model

$`\kappa _i^{}=p_i+{\displaystyle \frac{q_i}{N^2}}+{\displaystyle \frac{r_i}{N^4}},`$ (29)
with the parameters for each excited state given in Table 2.

These results should be compared again with those of Lucini and Teper who obtain large $`N`$ glueball masses in units of $`Ng^2/a`$ in the $`0^{}`$ sector with linear fits given by

$`0^{}:`$ $`1.176(14)+{\displaystyle \frac{0.14(20)}{N^2}}`$
$`0^{}:`$ $`1.535(28){\displaystyle \frac{0.35(35)}{N^2}}`$
$`0^{}:`$ $`1.77(13)+{\displaystyle \frac{0.24(161)}{N^2}}.`$ (30)
The result for $`0^{}`$ is from Ref. Teper, 1999. All fits were obtained using $`3N6`$ data. The $`N\mathrm{}`$ limit of our lowest lying massgap extracted from the large $`\beta `$ plateaux is consistent with Teper’s $`0^{}`$ result.

Having considered spin 0 states we now move on to spin 2, the only other spin accessible when using a basis of rectangular states. For the case of spin 2 the scaling behaviour of the massgaps is significantly worse than for spin 0. More troublesome is the fact that the convergence of the $`2^{++}`$ massgaps with increasing $`l_{\mathrm{max}}`$ is slower than the case of $`0^{}`$. The situation for the $`2^{}`$ state is markedly better with the convergence of massgaps with increasing $`l_{\mathrm{max}}`$ being no different to that of $`0^{++}`$. For this reason we concentrate on the $`2^{}`$ sector here.

The $`2^{}`$ states produce less precise scaling behaviour than the spin 0 states. However large $`\beta `$ plateaux appear for each of the five lowest lying massgaps. We use these regions to estimate their respective continuum limits. For these states small $`\beta `$ minima do not appear. The lowest lying $`2^{}`$ massgap is shown in Fig. 16 for a range of $`N`$. Once again we observe that in units of $`Ng^2/a`$ the massgaps do not depend strongly on $`N`$ and that in the scaling region they appear to approach a finite limit. Similar observations can be made for the four higher mass states obtained, although the scaling behaviour worsens as the mass of the state increases. The estimated continuum limit values extracted are plotted as functions of $`1/N^2`$ in Figs. 17 and 18. The dashed lines show fits to the model

$`\theta _i^{}=p_i+{\displaystyle \frac{q_i}{N^2}}+{\displaystyle \frac{r_i}{N^4}}.`$ (31)
The best fit parameters for each excited state are given in Table 3.

Once again these results should be compared with those of Lucini and Teper who obtain large $`N`$ glueball masses in units of $`Ng^2/a`$ in the $`2^{}`$ sector with linear fits derived from $`3N6`$ data given by

$`2^{}:`$ $`1.615(33){\displaystyle \frac{0.10(42)}{N^2}}`$
$`2^{}:`$ $`1.87(12){\displaystyle \frac{0.37(200)}{N^2}}.`$ (32)
The $`2^{}`$ result is from Ref. Teper, 1999. The $`N\mathrm{}`$ results presented here are significantly higher. This time no correspondence between our results and those of Lucini and Teper can be obtained.

We finish this section by presenting the mass spectra obtained for the $`0^{++}`$, $`0^{}`$ and $`2^{}`$ sectors. The results are summarised by the plots in Figs. 1921.

### III.5 An empirical observation

In this subsection we present an empirical observation that could possibly be useful in the construction of simple glueball models. Since Lucini and Teper include many small area Wilson loops in their construction of excited states, we assume that they calculate the lowest few glueball states without omission. We assume the calculations presented here give higher mass glueball states. Presumably, the correct enumeration of these states is not possible without the inclusion of nonrectangular states in the minimisation basis. For example, the third lowest mass state in a given sector calculated here could be the true seventh lowest mass state. However, when our results are combined with those of Lucini and Teper an interesting empirical observation can be made. We can choose a labelling of the $`0^{++}`$, $`0^{}`$ and $`2^{}`$ excited states such that the large $`N`$ limit of their masses lie on a straight line as shown in Fig. 22. We have included the results of Lucini and Teper for $`0^{}`$, $`0^{}`$, $`2^{}`$ and $`2^{}`$ which do not correspond to any of the states calculated here. The straight lines are fits to the model,

$`m_n(J^{PC})=\gamma _1(2n+\gamma _2)`$ (33)
where $`\gamma _1`$ and $`\gamma _2`$ are parameters and $`\gamma _2`$ is restricted to integer values. For the $`J^{PC}`$ states considered we obtain the following best fit models

$`m_n(0^{++})`$ $`=`$ $`(0.256\pm 0.002)(2n+1)`$
$`m_n(0^{})`$ $`=`$ $`(0.151\pm 0.002)(2n+5)`$
$`m_n(2^{})`$ $`=`$ $`(0.1495\pm 0.0008)(2n+7).`$ (34)
Here $`n1`$ labels the $`n`$-th lowest mass state. We note that the fit improves in accuracy as $`n`$ is increased. It is interesting to note that the slopes of the $`0^{}`$ and $`2^{}`$ models are consistent suggesting that the constant of proportionality does not depend on $`J`$. Another interesting observation can be made by recasting the models in the form

$`m_n(0^{++})`$ $`=`$ $`(0.256\pm 0.002)(2n+0+1)`$
$`m_n(0^{})`$ $`=`$ $`(0.151\pm 0.002)(2n+4+1)`$
$`m_n(2^{})`$ $`=`$ $`(0.1495\pm 0.0008)(2n+6+1).`$ (35)
We notice a similarity with the two dimensional harmonic oscillator spectrum,

$`E_n2n+J+1,`$ (36)
when we take into consideration the ambiguity modulo 4 of spin identification on the lattice. With this in mind we propose a simple model for the $`J^{PC}`$ spectrum

$`m_n(J^{PC})=\zeta _{PC}\left[2n+\gamma (J^{PC})+1\right],`$ (37)
where $`\zeta _{PC}`$ is a spin independent parameter and $`\gamma (J^{PC})`$ is an integer for which $`\gamma (J^{PC})=J\mathrm{mod}\mathrm{\hspace{0.17em}4}`$. From Eq. (34) we have $`\zeta _{}0.15`$ and $`\zeta _{++}=0.256\pm 0.002`$. To check this simple model we can attempt to predict the lowest lying states obtained by Lucini and Teper in the sectors that have not been considered in this paper. We start with the $`2^{++}`$ sector in which Teper and Lucini obtain $`1.359(12)`$ and $`1.822(62)`$ for the $`N\mathrm{}`$ limit of the $`2^{++}`$ and $`2^{++}`$ glueball masses in units of $`Ng^2/a`$. The predictions of Eq. (37), with $`\gamma (2^{++})=2`$ and $`\zeta _{++}=0.256`$, are 1.28 and 1.79 for the 1st and 2nd lowest glueball masses in units of $`Ng^2/a`$. As we would expect the prediction is better for the higher mass state. We can also consider the $`1^{++}`$ sector which has an equivalent spectrum to $`1^+`$ due to the phenomenon of parity doubling Teper (1999). Lucini and Teper obtain $`1.98(8)`$ for the mass of the $`1^{++}`$ glueball in units of $`Ng^2/a`$. The model of Eq. (37), with $`\gamma (1^{++})=5`$ and $`\zeta _{++}=0.256`$ gives 2.048. In the $`1^{}`$ sector the agreement is not as good, with Lucini and Teper obtaining $`1.85(13)`$ for the mass of the $`1^{}`$ glueball and the model of Eq. (37), with $`\gamma (1^{})=9`$ and $`\zeta _{}=0.15`$, giving 1.812.

It will be interesting to see if the model suggested here stands up to further calculations with an extended minimisation basis and in other $`J^{PC}`$ sectors. The analytic techniques presented here are at a great advantage to the standard Monte Carlo techniques of Lagrangian LGT for the purpose of testing glueball models. High order excited states are readily accessible; in a variational calculation with $`s`$ states in the minimisation basis, $`s`$ glueball states are accessible. This is in contrast to competing Lagrangian calculations in which only $`3`$ states are currently accessible in some $`J^{PC}`$ sectors.

### III.6 Discussion

In this section we have calculated variational mass spectra for pure SU($`N`$) gauge theory with a simple basis of rectangular states in 2+1 dimensions. Such a basis is easy to work with computationally and is therefore a good starting point. However to accurately explore the pure gauge mass spectra additional states need to be included. Perhaps a more suitable basis would be the set of all states with an area less than some maximum value which would define the order of the calculation. The next stage beyond this would be to include higher representation states in the form of Wilson loops in which some links are covered more than once. In this way multiple glueball states could presumably be explored. The inclusion of additional small area states in the basis, in particular states which are not symmetric about reflections in coordinate axes, would also allow the calculation of spin $`1`$ glueball masses and also massgaps in the $`P=C`$ sector.

It is important to point out that the inclusion of additional states in our minimisation basis does not present a significant challenge. The only complication would be in counting the possible overlaps of particular states. While more difficult than the counting required here, the process could presumably be automated using symbolic programing and techniques from graph theory. The analytic techniques used here would still be applicable until higher representation states were included in the minimisation basis. At that stage further character integrals would be required to handle the additional integrals arising in the calculation.

A final extension of the method presented here would be to make use of an improved ground state; a ground state which includes additional Wilson loops in the exponent. However, the advantage of using analytic techniques would then be lost, unless new technology was developed for tackling the resulting integrals. One would need to calculate the required matrix elements via Monte Carlo simulation on small lattices, losing the advantage of working on an infinite lattice with analytic expressions that we have here.

## IV Conclusion

In this paper we have applied the analytic techniques developed in Ref. Carlsson and McKellar, 2003, in a study of the large $`N`$ glueball mass spectrum in 2+1 dimensions.

In Section III we calculated glueball masses at finite $`N`$, with $`N`$ as large as 25 in some cases. This allowed accurate $`N\mathrm{}`$ extrapolations to be made. Evidence of leading order $`1/N^2`$ finite $`N`$ corrections to the glueball masses was obtained, confirming a specific prediction of large $`N`$ gauge theory. The interpretation of the possible scaling regions was discussed and agreement with the Lagrangian study of Lucini and Teper was obtained in some cases in the $`N\mathrm{}`$ limit. The discrepancies are attributable to our use of only rectangular states in the minimisation basis. Without the inclusion of small area, nonrectangular states, it is likely that some low energy states are inaccessible. Further work will fill in the incomplete spectra presented here. Interesting empirical observations were made when the results presented here were combined with those of Lucini and Teper; the enumeration of excited states can be chosen so that the glueball mass spectrum has the structure of a two dimensional harmonic oscillator. To develop this observation into a model would require a more complete calculation of the mass spectrum with nonrectangular states in the minimisation basis and possibly a more complicated vacuum trial state.

Clearly much work remains to be done before an accuruate picture of the 2+1 dimensional pure SU($`N`$) gauge theory spectrum is achieved within a Hamiltonian variational approach. We have demonstrated however that the analytic techniques of Hamiltonian LGT can be used for massgap calculations with $`N`$ as large as 25 on a desktop computer. This is significantly closer to the $`N\mathrm{}`$ limit than is currently possible with the Lagrangian approach using supercomputers. The inclusion of additional states in the minimisation basis, while not presenting a major challenge, will allow a thorough study of the 2+1 dimensional pure SU($`N`$) mass spectrum at least up to $`N=25`$.

###### Acknowledgements.

We wish to acknowledge useful and interesting discussions with J. A. L. McIntosh and L. C. L. Hollenberg.


text:
# Low-frequency measurement of the tunneling amplitude in a flux qubit

On leave from \]Department of Solid State Physics, Comenius University, SK-84248 Bratislava, Slovakia. On leave from \]Inst. of Microelectronic Technology, Russian Academy of Science, 142432 Chernogolovka, Russia. On leave from \]Novosibirsk State Technical University, 20 K. Marx Ave., 630092 Novosibirsk, Russia. On leave from \]B. Verkin Inst. for Low Temperature Physics and Engineering, 310164 Kharkov, Ukraine.

## Abstract

We have observed signatures of resonant tunneling in an Al three-junction qubit, inductively coupled to a Nb $`LC`$ tank circuit. The resonant properties of the tank oscillator are sensitive to the effective susceptibility (or inductance) of the qubit, which changes drastically as its flux states pass through degeneracy. The tunneling amplitude is estimated from the data. We find good agreement with the theoretical predictions in the regime of their validity.

Several groups, using different devices, have by now established that superconductors can behave as macroscopic quantum objects.macro ; Wal00 ; Rabi  These are natural candidates for a qubit, the building block of a quantum computer.

Qubits are effectively two-level systems with time-dependent parameters. One of them is a superconducting loop with low inductance $`L`$, including three Josephson junctions (a 3JJ qubit).Mooij99  Its potential energy, $`U=_{j=1}^3E_{\mathrm{J}j}(\varphi _j)`$, depends on the Josephson phase differences $`\varphi _j`$ across the junctions. Due to flux quantization $`_{j=1}^3\varphi _j=2\pi \mathrm{\Phi }_\mathrm{x}/\mathrm{\Phi }_0`$ (with $`\mathrm{\Phi }_\mathrm{x}`$ the external magnetic flux and $`\mathrm{\Phi }_0=h/2e`$ the flux quantum), only two $`\varphi _j`$’s are independent.

For suitable parameters, $`U(\varphi _1,\varphi _2)`$ has two minima corresponding to qubit states $`\mathrm{\Psi }^\mathrm{l}`$ and $`\mathrm{\Psi }^\mathrm{r}`$, carrying opposite supercurrents around the loop. These become degenerate for $`\mathrm{\Phi }_\mathrm{x}=\frac{1}{2}\mathrm{\Phi }_0`$. The Coulomb energy $`E_C`$ ($`e^2/2C`$, with $`C`$ the capacitance of junction 1) introduces quantum uncertainty in the $`\varphi _j`$. Hence, near degeneracy the system can tunnel between the two potential minima. (Since $`E_CE_\mathrm{J}E_{\mathrm{J1}}`$, we deal with a *flux* qubit; $`E_CE_\mathrm{J}`$ yields a *charge* qubit. Coherent tunneling was demonstrated in both.)

In the basis $`\{\mathrm{\Psi }^\mathrm{l},\mathrm{\Psi }^\mathrm{r}\}`$ and near $`\mathrm{\Phi }_\mathrm{x}=\frac{1}{2}\mathrm{\Phi }_0`$, the qubit can be described by the Hamiltonian

$$H(t)=ϵ(t)\sigma _z\mathrm{\Delta }\sigma _x;$$
(1)
$`\mathrm{\Delta }`$ is the tunneling amplitude. At bias $`ϵ=0`$ the two lowest energy levels of the qubit anticross \[Fig. 1(a)\], with a gap of $`2\mathrm{\Delta }`$. Increasing $`ϵ`$ slowly enough, the qubit can adiabatically transform from $`\mathrm{\Psi }^\mathrm{l}`$ to $`\mathrm{\Psi }^\mathrm{r}`$, staying in the ground state $`E_{}`$. Since $`dE_{}/d\mathrm{\Phi }_\mathrm{x}`$ is the persistent loop current, the curvature $`d^2E_{}/d\mathrm{\Phi }_\mathrm{x}^2`$ is related to the qubit’s susceptibility. Hence, near degeneracy the latter will have a peak, with a width given by $`|ϵ|\mathrm{\Delta }`$.Greenberg02b  We present data demonstrating such behavior in an Al 3JJ qubit.

Our technique is similar to rf-SQUID readout.Silver67 ; Ilichev01  The qubit loop is inductively coupled to a parallel resonant tank circuit \[Fig. 1(b)\]. The tank is fed a monochromatic rf signal at its resonant frequency $`\omega _\mathrm{T}`$. Then both amplitude $`v`$ and phase shift $`\chi `$ (with respect to the bias current $`I_\mathrm{b}`$) of the tank voltage will strongly depend on (A) the shift in resonant frequency due to the change of the effective qubit inductance by the tank flux, and (B) losses caused by field-induced transitions between the two qubit states. Thus, the tank both applies the probing field to the qubit, and detects its response.

The output signal depends on the tank’s quality factor $`Q`$. Using superconducting coil, values as high as $`Q10^3`$ can be obtained, leading to high readout sensitivity, e.g., in rf-SQUID magnetometers.Danilov  Such a tank can therefore be used to probe phase qubits.Greenberg02a  For small $`L`$, the results are summarized byGreenberg02b

$$v=I_0\omega _\mathrm{T}L_\mathrm{T}Q/\sqrt{1+(2Q\xi )^2},$$
(2)
$$\mathrm{tan}\chi =2Q\xi ,$$
(3)
$$\xi (v,f_\mathrm{x})=\frac{k^2L}{2\mathrm{\Phi }_0^2}_0^{2\pi }\frac{d\varphi }{\pi }\mathrm{cos}^2\varphi \frac{d^2E_{}(f)}{df^2},$$
(4)
$$f=f_\mathrm{x}+\frac{Mv}{\omega _\mathrm{T}L_\mathrm{T}\mathrm{\Phi }_0}\mathrm{sin}\varphi ,$$
(5)
where $`f_\mathrm{x}=\mathrm{\Phi }_\mathrm{x}/\mathrm{\Phi }_0\frac{1}{2}`$, $`I_0`$ is the bias-current amplitude, and $`k=M/\sqrt{LL_\mathrm{T}}`$ is the tank–qubit coupling coefficient, with $`M`$ ($`L_\mathrm{T}`$) the mutual (tank) inductance. The ground-state curvature istypo

$$\frac{d^2E_{}}{df^2}=\frac{E_\mathrm{J}^2\mathrm{\Delta }^2\lambda ^2}{(E_\mathrm{J}^2\lambda ^2f^2+\mathrm{\Delta }^2)^{3/2}},$$
(6)
where $`\lambda (\alpha ,g)`$ (with $`g=E_\mathrm{J}/E_C`$) is the conversion factor in $`ϵ=E_\mathrm{J}\lambda f`$.lambda  If $`I_0`$ vanishes, $`\xi =\frac{1}{2}k^2Ld^2E_{}/d\mathrm{\Phi }_\mathrm{x}^2`$ becomes an external parameter accounting for the qubit susceptibility coupled to the tank. For finite $`I_0`$, this has to be averaged over a bias cycle $`0<\varphi <2\pi `$. The resulting integral (4) turns out to involve a weight $`\mathrm{cos}^2\varphi `$, since the effective time-dependent coupling is $`(k\dot{f})^2`$ \[cf. the $`\varphi `$-derivative of Eq. (5)\], proportional to the square of the voltage the tank induces in the qubit. The resulting equations are coupled and nonlinear, but readily solved numerically.

For the tank, we prepared a square-shape Nb pancake coil on an oxidized Si substrate. The line width of the 20 windings was 2 $`\mu `$m, with a 2 $`\mu `$m spacing. Predefined alignment marks allow placing a qubit in the center. For flexibility, only the coil was made lithographically; an external capacitance $`C_\mathrm{T}`$ is used to change $`\omega _\mathrm{T}`$ in the range 5–35 MHz. For the selected tank ($`L_\mathrm{T}50`$ nH, $`C_\mathrm{T}470`$ pF), we obtained $`\omega _\mathrm{T}/2\pi =32.675`$ MHz and $`Q725`$ from the voltage–frequency characteristic.

The 3JJ qubit structure was manufactured out of Al by conventional shadow evaporation. The area of two of the junctions was estimated using electron microscopy as $`190\times 650`$ nm2 while one is smaller, so that $`\alpha E_{\mathrm{J3}}/E_{\mathrm{J1},2}0.8`$. The critical current was determined by measuring an rf-SQUID prepared on the same chipIlichev01  as $`I_\mathrm{c}=2eE_\mathrm{J}/\mathrm{}380`$ nA. With $`E_C/h3`$ GHz, one finds $`g60`$ and $`\lambda 4.4`$. The loop area was 90 $`\mu \mathrm{m}^2`$, with $`L=38`$ pH. We measured $`v`$ by a three-stage cryogenic amplifier, placed at $`2`$ K and based on commercial pseudomorphic high electron mobility transistors. It was slightly modified from the version in Ref. Oukhansky02,  to decrease its back-action on the qubit. The input-voltage noise was $`<0.6`$ nV$`/\sqrt{\mathrm{Hz}}`$ in the range 1–35 MHz. The noise temperature was $`300`$ mK at 32 MHz. The effective qubit temperature due to the amplifier’s back-action should be considerably lower because of the small $`k210^2`$.

The $`\chi (f_\mathrm{x})`$ curves measured at various $`I_0`$ and a mixing-chamber temperature $`T=10`$ mK are shown in Fig. 2. The narrow dip at $`f_\mathrm{x}=0`$ directly corresponds to the one in Eq. (6), in line with the qualitative picture below Eq. (1). With device parameters as above, all quantities in Eqs. (2)–(6) are known, but $`\mathrm{\Delta }`$ only in principle: its exponential sensitivity to $`\alpha `$ and especially $`g`$ makes it notoriously hard to calculate *a priori*. Hence, it is treated as a free parameter; calculated curves for the best fit $`\mathrm{\Delta }/h=650`$ MHz are shown in the inset. For the largest $`I_0`$ the experimental and theoretical curves disagree, for the rapid change of $`\mathrm{\Phi }_\mathrm{x}`$ then leads to Landau–Zener transitionsLZ ; andrei  suppressing the dip.

The $`T`$-dependence of $`\chi `$ is shown in Fig. 3. For increasing $`T`$ the dip’s amplitude decreases while, strikingly, its width is unchanged \[Fig. 3(c)\]. Both are a simple manifestation of the Hamiltonian (1) yielding $`\sigma _z=(ϵ/\mathrm{\Omega })\mathrm{tanh}(\mathrm{\Omega }/k_\mathrm{B}T)`$,adia  $`\mathrm{\Omega }=\sqrt{ϵ^2+\mathrm{\Delta }^2}`$. This result of equilibrium statistics of course assumes that the $`t`$-dependence of $`ϵ(t)`$ is adiabatic. However, it *does* remain valid if the full (Liouville) evolution operator of the qubit would contain standard Bloch-type relaxation and dephasing terms (which indeed are not probedGreenberg02b ) in addition to the Hamiltonian dynamics (1), since the fluctuation–dissipation theorem guarantees that such terms do not affect equilibrium properties. Normalized dip amplitudes are shown vs $`T`$ in Fig. 3(b) together with $`\mathrm{tanh}(\mathrm{\Delta }/k_\mathrm{B}T)`$, for $`\mathrm{\Delta }/h=650`$ MHz *independently* obtained above from the low-$`T`$ width. The good agreement strongly supports our interpretation, and is consistent with $`\mathrm{\Delta }`$ being $`T`$-independent in the relevant range.Han02  Of course, for higher $`T`$ the dip will wash out; we observe a width $`T^3`$ above a crossover temperature $`225`$ mK. For $`T`$ of this order, deviations from the two-state model can be expected, especially for $`f_\mathrm{x}0`$. This behavior outside the qubit regime has not been pursued.

In conclusion, we have observed resonant tunneling in a macroscopic superconducting system, containing an Al flux qubit and a Nb tank circuit. The latter played dual control and readout roles. The impedance readout technique allows direct characterization of some of the qubit’s quantum properties, *without* using spectroscopy.Wal00 ; Rabi  In a range 50$``$200 mK, the *effective* qubit temperature has been verified \[Fig. 3(b)\] to be the same as the mixing chamber’s (after $`\mathrm{\Delta }`$ has been determined at low $`T`$), which is often difficult to confirm independently.

MHSA and AMvdB are grateful to A.Yu. Smirnov and A.M. Zagoskin for fruitful discussions, and to P.C.E. Stamp for the remark on effective thermometry. MG wants to acknowledge partial support by the Slovak Grant Agency VEGA (Grant No. 1/9177/02).