text:
# Josephson effects in dilute Bose-Einstein condensates

## Abstract

We propose an experiment that would demonstrate the “dc” and “ac” Josephson effects in two weakly linked Bose-Einstein condensates. We consider a time-dependent barrier, moving adiabatically across the trapping potential. The phase dynamics are governed by a “driven-pendulum” equation, as in current-driven superconducting Josephson junctions. At a critical velocity of the barrier (proportional to the critical tunneling current), there is a sharp transition between the “dc” and “ac” regimes. The signature is a sudden jump of a large fraction of the relative condensate population. Analytical predictions are compared with a full numerical solution of the time dependent Gross-Pitaevskii equation, in an experimentally realistic situation.

The Josephson effects (JE’s) are a paradigm of the phase coherence manifestation in a macroscopic quantum system . Observed early on in superconductors , JE’s have been demonstrated in two weakly linked superfluid 3He-B reservoirs . Weakly interacting Bose-Einstein condensate (BEC) gases  provide a further (and different) context for JE’s. Indeed, magnetic and optical traps can be tailored and biased (by time-dependent external probes) with high accuracy , allowing the investigation of dynamical regimes that might not be accessible with other superconducting/superfluid systems. The macroscopic BEC’s coherence has been demonstrated by interference experiments , and the first evidence of coherent tunneling in an atomic array, related to the “ac” JE, has been recently reported .

A superconducting Josephson junction (SJJ) is usually biased by an external circuit that typically includes a current drive $`I_{ext}`$. The striking signatures of the Josephson effects in SJJ are contained in the voltage-current characteristic ($`V`$-$`I_{ext}`$), where usually one can distinguish between the superconductive branch or “dc”-branch ($`V=0`$, $`I_{ext}0`$), and the resistive branch or “ac”-branch ($`VRI_{ext}`$) (see for example ). External circuits and current sources are absent in two weakly linked Bose condensates and the Josephson effects have been related, so far, with coherent density oscillations between condensates in two traps or between condensates in two different hyperfine levels . This collective dynamical behavior is described by a non-rigid pendulum equation , predicting a new class of phenomena not observable with SJJ’s.

Now the following question arises: can two weakly linked condensates exhibit the analog of the voltage-current characteristic in SJJ? Although BECs are obviously neutral, the answer is positive. A dc current-biased SJJ can be simulated by considering a tunneling barrier moving with constant velocity across the trap. At a critical velocity of the barrier a sharp transition between the “dc” and “ac” (boson) Josephson regimes occurs. This transition is associated with a macroscopic jump in the population difference, that can be easily monitored experimentally by destructive or non-destructive techniques.

In the following we will briefly introduce the phenomenological equations of the resistively shunted junction (RSJ) model for the SJJ. We will describe the corresponding experiment for two weakly linked BECs and show that the relevant equations are formally equivalent to the RSJ equations. Then we compare the analytical results with a numerical integration of the Gross-Pitaevskii equation in a realistic 3D setup.

In the RSJ model, SJJ is described by an equivalent circuit  in which the current balance equation is

$$I_c\mathrm{sin}(\theta )+GV+C\dot{V}=I_{ext}$$
(1)
where $`I_c`$ is the upper bound of the Josephson supercurrent $`I`$ (which is represented, in the ideal case, by the sinusoidal current-phase relation $`I=I_c\mathrm{sin}(\theta )`$); $`G`$ is an effective conductance (offered by the quasiparticles and the circuit shunt resistor), and $`C`$ is the junction capacitance. The voltage difference $`V`$ across the junction is related to the relative phase $`\theta `$ by

$$\dot{\theta }=2eV/\mathrm{}.$$
(2)
In the low conductance limit $`G\omega _pC`$ where $`\omega _p=\sqrt{2eI_c/\mathrm{}C}`$ is the Josephson plasma frequency, combining equations (1) and (2) leads to the “driven pendulum” equation

$$\ddot{\theta }=\omega _p^2\frac{}{\theta }U\left(\theta \right)$$
(3)
where $`U`$ is the tilted “washboard” potential:

$$U\left(\theta \right)=1\mathrm{cos}(\theta )+i\theta $$
(4)
with $`i=I_{ext}/I_c`$. This equation describes the transient behavior before the stationary dissipative behavior is reached (resistive branch). If we start from equilibrium, with $`i=0`$, and increase adiabatically the current, no voltage drop develops until the critical value $`i=1`$ is reached (neglecting secondary quantum effects). At this point $`V`$ continuously develops until a stationary asymptotic dissipative behavior is reached in a time scale approximately of order $`C/G`$. Similar phenomenology may occur in BECs and we will derive equations formally identical to Equations (3) and (4).

A weak link between two condensates can be created by focusing a blue-detuned far-off-resonant laser sheet into the center of the magnetic trap . The weak link can be tailored by tuning the width and/or the height of the laser sheet. Raman transitions between two condensates in different hyperfine levels provide a different weak link , in analogy with the “internal Josephson effect” observed in $`70`$s with $`{}_{}{}^{3}HeA`$ .

Here we consider a double well potential in which the laser sheet slowly moves across the magnetic trap with velocity $`v`$ (but our framework can be easily adapted to investigate the internal Josephson effect). In the limit of very low $`v`$, the two condensates remain in equilibrium, i.e. in their instantaneous ground state, because of the non-zero tunneling current that can be supported by the barrier. In fact, an average net current, proportional to the velocity of the laser sheet, flows through the barrier, sustained by a constant relative phase between the two condensates. This keeps the chemical potential difference between the two subsystems locked to zero, as in the SJJ dc-branch. However, the superfluid component of the current flowing through the barrier is bounded by a critical value $`I_c`$. As a consequence there exists a critical velocity $`v_c`$, above which a non-zero chemical potential difference develops across the junction. This regime is characterized by a running-phase mode, and provides the analog of the ac-branch in SJJ’s.

The ”dc” and ”ac” BEC regimes are governed by a phase-equation similar to the current-driven pendulum equations (3) and (4). Such equations together with the sinusoidal current-phase relation $`I=I_c\mathrm{sin}(\theta )`$ describe the phase difference and current dynamics. The dimensionless current $`i`$ is related to the barrier velocity by

$$i=v/v_c$$
(5)
with the critical velocity $`v_c`$ given by

$$v_c=\frac{\mathrm{}\omega _p^2}{F}$$
(6)
where $`F`$ is to a good approximation represented by double the average force exerted by the magnetic trap on single atoms in one well.

Equations (3)-(6) can be derived by a time-dependent variational approximation and have also been verified, as we discuss below, by the full numerical integration  of the Gross-Pitaevskii equation (GPE) . The GPE describes the collective dynamics of a dilute Bose gas at zero temperature:

$$i\mathrm{}\frac{}{t}\mathrm{\Psi }=\left[H_0\left(t\right)+g|\mathrm{\Psi }|^2\right]\mathrm{\Psi }$$
(7)
where $`H_0\left(t\right)=\frac{\mathrm{}^2}{2m}^2+V_{ext}(𝐫,t)`$ is the non interacting Hamiltonian and where $`g=4\pi \mathrm{}^2a/m`$, with $`a`$ the scattering length and $`m`$ the atomic mass. The order parameter $`\mathrm{\Psi }=\mathrm{\Psi }(𝐫,t)`$ is normalized as $`𝑑𝐫|\mathrm{\Psi }(𝐫,t)|^2=N`$, with $`N`$ the total number of atoms. The external potential is given by the magnetic trap and the laser barrier $`V_{ext}(𝐫,t)=V_{trap}\left(𝐫\right)+V_{laser}(z,t)`$. We consider a harmonic, cylindrically symmetric trap $`V_{trap}\left(𝐫\right)=\frac{1}{2}m\omega _r^2\left(x^2+y^2\right)+\frac{1}{2}m\omega _0^2z^2`$ where $`\omega _r`$ and $`\omega _0`$ are the radial and longitudinal frequency, respectively. The barrier is provided by a Gaussian shaped laser sheet, focused near the center of the trap $`V_{laser}\left(z\right)=V_0\mathrm{exp}\left((zl_z)^2/\lambda ^2\right)`$ with the coordinate $`l_z(t)`$ describing the laser motion and $`v=dl_z/dt`$ its velocity.

The equations (3) to (6) can be derived by solving variationally the GPE using the ansatz: $`\mathrm{\Psi }(𝐫,t)=c_1(t)\psi _1\left(𝐫\right)+c_2(t)\psi _2\left(𝐫\right)`$, where $`c_n=\sqrt{N_n(t)}\mathrm{exp}\left(i\theta _n(t)\right)`$ are complex time-dependent amplitudes of the left $`n=1`$ and right $`n=2`$ condensates (see also ). The trial wave functions $`\psi _{1,2}\left(𝐫\right)`$ are orthonormal and can be interpreted as approximate ground state solutions of the GPE of the left and right wells. The equations of motion for the relative population $`\eta =(N_2N_1)/N`$ and phase $`\theta =\theta _2\theta _1`$ between the two symmetric traps are

$`\mathrm{}\dot{\eta }`$ $`=`$ $`(2E_J/N)\sqrt{1\eta ^2}\mathrm{sin}\left(\theta \right),`$ (8)
$`\mathrm{}\dot{\theta }`$ $`=`$ $`Fl_z(t){\displaystyle \frac{2E_J}{N}}{\displaystyle \frac{\eta }{\sqrt{1\eta ^2}}}\mathrm{cos}\left(\theta \right){\displaystyle \frac{NE_c}{2}}\eta ,`$ (9)
where $`E_c=2g𝑑𝐫\psi _1(𝐫)^4`$ is the variational analog of the capacitive energy in SJJ, while $`E_J=N𝑑𝐫\psi _1(𝐫)\left[H_0+gN\psi _1^2(𝐫)\right]\psi _2(𝐫)`$ is the Josephson coupling energy. The current-phase relation $`I=I_c\sqrt{1\eta ^2}\mathrm{sin}(\theta )`$ is directly related to Eq. (8) where the critical current is given by $`I_c=E_J/\mathrm{}`$. $`Fl_z(t)`$ represents the contribution to the chemical potential difference in the two wells due to the laser displacement $`l_z`$ (after linearizing in $`l_z`$), and where $`F=𝑑𝐫\left(\psi _1(𝐫)^2\psi _2(𝐫)^2\right)\frac{}{l_z}V_{laser}m\omega _0^2𝑑𝐫z\left(\psi _1(𝐫)^2\psi _2(𝐫)^2\right)`$. The above variational method provides a simple and useful interpolating scheme between the low interacting limit $`N^2E_cE_J`$ and the opposite limit $`N^2E_cE_J`$. In the last case, and with $`\eta 1`$, we recover the driven-pendulum phase equation (3) and the critical velocity relations (5) and (6) with $`\mathrm{}\omega _p=\sqrt{E_JE_c}`$. In particular, it is legitimate to consider the Josephson coupling as a perturbation, with the the phase dynamics entirely determined by the difference in the chemical potentials $`\mu _1(N_1,l_z)`$ and $`\mu _2(N_2,l_z)`$ in the two wells. In this case $`E_c`$ corresponds to $`2\left(\mu _1/N_1\right)_{l_z}`$ and $`\mathrm{}^2\omega _p^2=E_J\left(\mu _1/N_1\right)_{l_z}`$. The critical velocity is proportional to the critical current: $`v_c=\left(\frac{dN_1}{dl_z}\right)^1I_c`$, with

$$\left(\frac{dN_1}{dl_z}\right)^1=\left(\frac{\mu _1}{l_z}\right)_{N_1}^1\left(\frac{\mu _1}{N_1}\right)_{l_z}$$
(10)
and $`\left(\mu _1/l_z\right)_{N_1}`$ being $`F/2`$ in Eq.(6). These derivatives can be computed numerically. In the Thomas-Fermi (TF) limit they reduce to

$$\left(\frac{\mu _1}{N_1}\right)_{l_z}=\frac{g}{V_{TF}}$$
(11)
and

$$\left(\frac{\mu _1}{l_z}\right)_{N_1}=\frac{1}{V_{TF}}_{V_{TF}}𝑑𝐫\frac{}{l_z}V_{laser}$$
(12)
where $`V_{TF}`$ is the volume of the region in which $`\mathrm{\Psi }_1`$ is different from zero (in the TF approximation).

We make the comparison of Eqs. (8) and (9) with a full numerical integration of the GPE in an experimentally realistic geometry relative to the limit $`N^2E_cE_J`$. In particular, we show that Eq. (6), derived in the limit of $`\eta 1`$, still remains a good approximation even for $`\eta 0.4`$. The details of the numerical calculation are given elsewhere .

We have considered the JILA setup, with $`N=5\times 10^4`$ Rb atoms in a cylindrically symmetric harmonic trap, having the longitudinal frequency $`\omega _0=50`$ s-1 and the radial frequency $`\omega _r=17.68`$ s-1. The value of the scattering length considered is $`a=58.19`$ $`\dot{A}`$. A Gaussian shaped laser sheet is focused in the center of the trap, cutting it into two parts. We assume that the (longitudinal) $`1/e^2`$ half-width of the laser barrier is $`3.5`$ $`\mu `$m and the barrier height $`V_0/\mathrm{}=650`$ s-1.

Although the lifetime of a trapped condensate can be as long as minutes, we have made a quite conservative choice, by considering a time scale on the order of one second. The possibility to perform experiments on a longer time-scale will improve the observability of the phenomena we are discussing. With this choice of time scale, that corresponds only to few plasma oscillations, an adiabatic increase of the velocity is not possible, therefore we proceed as follows. For $`t<0`$ the laser is at rest in the middle of the trap, $`l_z=0`$, and the two condensates are in equilibrium. For $`t>0`$ the laser moves across the trap, with constant velocity, and the relative atomic population is observed at $`t_f=1s`$. With this initial condition, which introduces small plasma oscillations in the relative population, it is expected, in absence of dissipation, to slightly reduce the critical current by the numerical factor $`0.725`$ (see the general properties of the driven pendulum equation ).

In Fig.1 we show the relative condensate population $`\eta =(N_2N_1)/N`$, calculated after $`1`$ second, for different values of the laser velocity $`v`$. The crosses are the results obtained with the full numerical integration of the time-dependent GPE (7). The dot-dashed line shows the equilibrium values $`\eta _{eq}`$ of the relative population calculated with the stationary GPE and with the laser at rest in the ”final” position $`l_z=vt_f`$. The displacement of $`\eta (t_f)`$ from $`\eta _{eq}`$ is a measure of the chemical potential difference, being $`\mathrm{\Delta }\mu =\mu _2\mu _1NE_c(\eta (t_f)\eta _{eq})/2`$.

For $`v<0.42\mu m/s`$, the atoms tunnel through the barrier in order to keep the chemical potential difference $`\mathrm{\Delta }\mu `$ locked around zero. The dc component of the tunneling current is accounted for by an averaged constant phase difference between the two condensates. This is the close analog of the dc Josephson effect in superconducting Josephson junctions. The small deviations between the dashed line and the crosses are due to the presence of plasma oscillations (induced by our initial condition). At $`v0.42\mu m/s`$ there is a sharp transition, connected with the crossover from the dc-branch to the ac-branch in SJJ. For $`v>0.42\mu m/s`$, the phase difference starts running and the population difference, after a transient time, remains on average fixed. A macroscopic chemical potential difference is established across the junction. In this regime ac oscillations in the population difference are observed. The frequency of such oscillations are approximatively given by $`\mathrm{\Delta }\mu (t)/\mathrm{}`$ (not visible in the figure).

The solid line of Fig.1 corresponds to the solutions of Eqs. (8) and (9) in which the value of the energy integrals $`E_cN/\mathrm{}=2.46ms^1`$ and $`E_J/N\mathrm{}=2.41\times 10^4ms^1`$ are chosen in order to give the correct value of $`\omega _p=2.44\times 10^2ms^1`$ and $`I_c=12.1ms^1`$. The values $`\omega _p`$, $`I_c`$ are calculated numerically studying the frequency of small oscillations around equilibrium and the current-phase relation, respectively. The force integral is $`F/\mathrm{}=1.060ms^1\mu m^1`$. The parameters $`\omega _p`$, $`I_c`$ and $`F`$ are calculated with the laser at rest ($`v=0`$) in $`l_z=0`$. Using these values in Eq. (6) and taking into account the reducing factor $`0.725`$ we obtain the value $`0.407`$ $`\mu ms^1`$ for the critical velocity, in agreement with the value observed in the simulation.

Small deviations between the variational solutions (full line in Fig.1) and the numerical results (crosses in Fig.1), above the critical velocity, are due to “level-crossing” effects. Numerical results  show that when the condensate ground state of the “upper” well is aligned with the excited collective dipole state in the “lower” well, a finite number of atoms go from the “upper” well to the “lower” well. Close to this tunneling resonance it is possible to control, by manipulating the barrier velocity below a fraction of $`v_c`$, the dc flux of atoms from the ground state condensate in the “upper” well to the longitudinal intrawell collective dipole mode of the condensate in the “lower” well. This effect is directly observable in the macroscopic longitudinal oscillations of the two condensates (at frequencies $`\omega _0`$).

Concerning a possible realization of the phenomenon described in this work, we note that for small barrier velocities $`v`$, the motion of the laser sheet with respect to the magnetic trap with velocity $`v`$ or, $`viceversa`$, the motion of the magnetic trap with velocity $`v`$, are equivalent, there being negligible corrections due to different initial accelerations.

Thus far we have discussed the zero temperature limit. At finite temperature dissipation can arise due to incoherent exchange of thermal atoms between the two wells. This can be described phenomenologically by including a term $`E_cG\dot{\theta }/\omega _p^2`$ in Eq. (3) where $`G`$ is the conductance. Dissipation will be negligible as long as the characteristic time scale $`(E_cG)^1(20G/\mathrm{})s`$ is bigger than the time scale of the experiment ($`1s`$).

To conclude we note that while it could be difficult to measure directly the plasma oscillations, since their amplitude is limited by $`\mathrm{\Delta }\eta <\frac{4}{N}\sqrt{\frac{E_J}{E_c}}`$, the macroscopic change in the population difference may be easily detected with standard techniques. Moreover the framework that we have discussed can be easily adapted to investigate the internal Josephson effect.

Our phenomenological equations are similar to the driven pendulum equation governing the Josephson effects in SJJs. As a consequence, within this framework we can study the “secondary quantum phenomena”, such as the Macroscopic Quantum Tunneling between different local minima of the washboard potential (see for instance ).

It is a pleasure to thank L. P. Pitaevskii, S. Raghavan and S. R. Shenoy for many fruitful discussions.


format: text

text:
# NICMOS Narrow-band Infrared Photometry of TW Hya Association Stars

## 1 Introduction

For more than a decade, the young star TW Hya has been an enigma since it lies in a region of sky apparently devoid of the raw materials to form stars, nearly 13 from the nearest dark cloud, yet it is unambiguously a classical T Tauri star (Rucinski & Krautter (1983)) surrounded by a great deal of cold dust (Weintraub, Sandell & Duncan (1989)) and gas (Zuckerman et al. (1995); Kastner et al. (1997)). Recently, the TW Hya mystery was solved: TW Hya, along with other T Tauri stars found in an area of $``$100 square degrees of the southern sky (de la Reza et al. (1989); Gregorio-Hetem et al. (1992)), compose a uniquely close association of young stars known as the TW Hya Association (Kastner et al. (1997)).

At a mean distance of only $``$55 pc, the TW Hya Association (hereafter TWA) is almost three times closer than the next nearest known region of recent star formation. Given the likely age ($``$10 MY) of the TWA, these stars could harbor very young planetary systems with fully formed giant planets or low mass, brown dwarf companions, and may still be surrounding by circumstellar disks. In fact, there is substantial evidence for circumstellar gas and dust around several of these stars (Weintraub, Sandell & Duncan (1989); Zuckerman & Becklin (1993); Zuckerman et al. (1995); Kastner et al. (1997)). The relative proximity and the absence of significant interstellar or intra-molecular cloud extinction in the direction of the TWA make the prospects for detecting substellar companions around these nearest T Tauri stars much better than the prospects for similar searches for young low mass companions around T Tauri stars in Taurus-Auriga, Chamaeleon, Lupus or Ophiuchus, the next closest regions of star formation. In a recent study of the TWA, Webb et al. (1999) identified a total of at least 17 sources as members of the TWA. In addition, Lowrance et al. (1999) and Webb et al. have reported the discovery of a likely low mass brown dwarf companion (M $``$ 0.02 $`M_{}`$) to TWA 5A (= CD$``$337795) in a combination of ground-based and Hubble Space Telescope (HST) observations. TWA 5B is found almost 2′′ from TWA 5A; thus, despite its relative physical proximity ($``$100 AU) to the primary, TWA 5B is amenable to spectroscopic and astrometric studies, uncontaminated by light from TWA 5A.

In this paper, we report results from imaging the fields around five stars in the TWA, including the TWA 5 system, using the Near Infrared Camera and Multi-Object Spectrometer (NICMOS) and the HST. The goal of this program was to search for companions around these stars. Our choice of three narrowband filters centered at 1.64, 1.90 and 2.15 $`\mu `$m was designed to enable us to identify cool and low surface gravity objects, including substellar mass companions, through their likely strong signatures of H2O absorption at 1.9 $`\mu `$m.

## 2 Observations

We obtained images of five star systems (Table 1) in the TWA using camera 1 (NIC1) and camera 2 (NIC2) of NICMOS between 1998 May 30 and July 12 (U.T.) Observations of each of the five stars were made identically. Using NIC1 and filter F164N, we imaged each target in a four position, spiral dither pattern, with an integration time per position of 33.894 s. Three images were obtained at each position for a total integration time of 406.73 s. We carried out identical observations using NIC1 and filter F190N, with an integration time per position of 43.864 s and a total integration time of 526.37 s. Switching to NIC2 and filter F215N, we again obtained three sets of four-position dithered image suites; however, for the F215N observations we changed the starting position for the dithered image suites in order to obtain a better median filtered image for subtraction of the thermal background. The integration time per position was 15.948 s for the NIC2 images, for a total integration time of 191.38 s.

## 3 Results

### 3.1 Imaging

We find no sources in any of our images other than the previously known five primaries and four secondaries, to limiting magnitudes of 18.3, 18.4 and 17.5 in the F164N, F190N and F215N images, respectively, at distances beyond $``$1$`^{\prime \prime }.`$3 at 1.64 and 1.90 $`\mu `$m and 2$`^{\prime \prime }.`$3 at 2.15 $`\mu `$m. The images of TWA 5 (Fig. 1) reveal how how easy it is to detect and image young, intermediate mass, brown dwarf companions around stars in the TWA, even without a coronagraph. In addition, all nine imaged objects appear as point sources (with FWHM of 0$`^{\prime \prime }.`$14, 0$`^{\prime \prime }.`$16, and 0$`^{\prime \prime }.`$18 in the F164N, F190N and F215N images, respectively), with no evidence (after deconvolutions performed with point spread functions \[PSFs\] generated using the software package Tiny Tim111http://scivax.stsci.edu/ krist/tinytim.html, direct subtractions of PSFs, and examinations of azimuthally averaged radial intensity profiles \[Fig. 2\]) of extended emission around any of them. Thus, although some of these stars appear to be surrounded by circumstellar material (e.g., TWA 1 is surrounded by a circumstellar disk of radius $``$3′′ that is viewed nearly face-on; Weinberger et al. (1999); Krist et al. (1999)), we conclude that these direct, narrow band images are insufficiently sensitive to image circumstellar disks around these stars.

### 3.2 Photometry

We report our photometry for these observations in Table 1. The factors used to convert from NICMOS count rates to absolute fluxes and magnitudes222http://www.stsci.edu/ftp/instrument\_news/NICMOS/NICMOS\_phot/keywords.html, version 1998, December 1 were 5.376665 $`\times `$ 10-5 Jy sec ADU-1 for the F164N filter, 4.866353 $`\times `$ 10-5 Jy sec ADU-1 for the F190N filter, and 3.974405 $`\times `$ 10-5 Jy sec ADU-1 for the F215N filter with zero point flux densities of 1033 Jy, 862 Jy, and 690 Jy, respectively. Photometry was obtained by measuring the total counts within a 0$`^{\prime \prime }.`$5 radius aperture and then applying a correction factor of 1.15 to compensate for the flux that falls outside of this radius333 http://www.stsci.edu/ftp/instrument\_news/NICMOS/nicmos\_doc\_phot.html.

Figure 3 shows how the photometry of TWA stars through the $`J`$, $`H`$, and $`K`$ broadband and the F164N, F190N and F215N narrowband filters relates to the near infrared spectral characteristics of late-type stars. Stars with lower effective temperatures have increasingly strong water absorption bands centered at $`1.4\mu `$m and $`1.9\mu `$m which are very effectively probed by this combination of broad and narrowband filters.

We find no systematic differences between the $`K`$ and F215N photometry (see Fig. 3, right panel), despite the slight difference in central wavelengths and large difference in bandwidth. On the other hand the stars are systematically brighter at F164N than at $`H`$-band, the most extreme case being TWA 5B. Finally, TWA 5B is the only star with an absolute flux that is clearly lower at 1.9 than at 2.15 $`\mu `$m.

The left hand panel of Fig. 3 clearly shows how water absorption in the 1.35-1.55 and 1.7 - 2.1 $`\mu `$m regions will strongly affect broad band $`H`$ measurements but will have no effect on observations with the F164N filter (see also the library of near-infrared spectra published by Lançon & Rocca-Volmerange 1992). Thus, the $`H`$-band and F164N observations reveal the presence of different amounts of water vapor absorption in the spectra of most of the stars in our sample.

### 3.3 Astrometry

For the four binary systems in our sample, we measured the intensity centroids for each binary component and transformed the cartesian positions on the array into offsets in Right Ascension and Declination of the secondaries from the primaries (Table 2). Except for TWA 8, for which the binary separation is such that the companion only appeared in the NIC2 images, the results presented in Table 2 are those obtained using only the NIC1 images since the spatial resolution is highest when using the NIC1 array. The NIC1-based results in Table 2 are the statistical average and standard deviations based on measurements of the F164N and F190N images. The pixel to RA and Dec conversions were done using the plate scale measurement ephemeris generated by the NICMOS instrument team444 http://www.stsci.edu/ftp/instrument\_news/NICMOS/nicmos\_doc\_platescale.html. Comparison of the results for the F164N and F190N images indicate that, in most cases, we can determine image separations to an accuracy of 0.02 pixels ($`<`$ 1 milli-arcsec).

Because our astrometric results are obtained from unocculted HST images and with the highest resolution camera in NICMOS, these results are much more precise than offsets previously reported for these binaries. They are, however, consistent with previous results (Table 2). In the case of TWA 5B, Lowrance, Weinberger & Schneider (1999) recently independently determined that the offsets reported in Lowrance et al. (1999) and Webb et al. (1999) have a sign error in RA; the corrected values are reported in Table 2.

## 4 TWA 5B

### 4.1 The age of the TWA

We have constructed an H-R diagram for the TWA (Fig. 4) using the pre-main sequence tracks of Baraffe et al. (1998). This H-R diagram is quite similar to that presented by Webb et al. (1999), which is based on the pre-main sequence tracks of D’Antona & Mazzitelli (1997); however, Fig. 4 appears to constrain the cluster age more tightly than does previous work on the TWA, presumably because of improved physics included in the Baraffe et al. tracks (see Baraffe et al. 1997 for a discussion). Specifically, virtually all of the stars, including TWA 5B, fall between the 3 and 10 MY isochrones. In comparison, the H-R diagrams of Webb et al. and Lowrance et al. (1999) indicate that the TWA stars have ages in the range 1–100 MY while Kastner et al. (1997) suggested that the likely age of the TWA stars is 10–30 MY, based on lithium studies (upper limit) and X-ray luminosities (lower limit). TWA 6 and TWA 9A, which lie together almost on the 30 MY isochrone, and TWA 9B, which falls near the 100 MY isochrone, are mild outliers in our and the Webb et al. HR diagrams and appear to be older than the other TWA stars (however, see Webb et al. for other possible explanations).

What other information do we have to constrain the ages of the TWA stars? Soderblom et al. (1998) used the lithium abundance to place an age range of 5–20 MY and a most probable age of 10 $`\pm `$ 3 MY on TWA 4 (HD 98800; EW(Li $`\lambda `$6708) = 0.36 Å); Stauffer, Hartmann, & Barrado y Navascues (1995) used the strength of the Li line to assign an upper limit of 9–11 MY to TWA 11B (HR 4796B) while Jayawardhana et al. (1998) assigned an isochronal age of 8 $`\pm `$ 3 MY to this star; and Webb et al. (1999) measured similar Li EW strengths for 14 of the 17 stars identified as members of the TWA and, on this basis, suggested that they are all less than $``$10 MY. The excellent agreement between the ages estimated from the Li EWs and those obtained from photometry and pre-main sequence evolutionary tracks suggests that the age of the TWA is well constrained to be in the range 5–15 MY.

### 4.2 Mass and evolutionary status of TWA 5A and 5B

TWA 5A is a M1.5 star (Webb et al. 1999) with $`T_{\mathrm{eff}}=3700\pm 150`$K (Leggett et al. 1996). The distance to the TWA 5 system is presently unknown but can be estimated as $`55\pm 9`$pc from the measured parallaxes of four members of the association (Webb et al. 1999). The range of distances is consistent with the approximate angular dimension of the association. With $`d=55\pm 9`$pc and $`K=6.8\pm 0.1`$, we find $`M_K=3.10\pm 0.41`$. By comparing these values of $`T_{\mathrm{eff}}`$ and $`M_K`$ with the evolution sequences of Baraffe et al. (1998), we find $`M=0.75\pm 0.15M_{}`$ and an age of 2.5 to 6 MY for TWA 5A, assuming it is a single, pre-main sequence star (Fig. 4). On the other hand, Webb et al. report that TWA 5A is suspected to be a spectroscopic binary. If we assume that TWA 5A is binary with equal mass components, the mass of each component decreases to $`0.7\pm 0.15M_{}`$ and the age range becomes 6 to 18 MY.

Since the Baraffe et al. (1998) sequence does not extend to substellar masses, we analyze the photometric measurements of TWA 5B with evolutionary models computed by Saumon & Burrows (unpublished). These models use the same interior physics as Saumon et al. (1996) and Burrows et al. (1997) with the distinction that the surface boundary condition is provided by the “NextGen” sequence of atmosphere models computed by Allard and Hauschildt for cool stars (Allard et al. 1996, Hauschildt, Allard & Baron 1999). The atmospheric structures provide a surface boundary condition for the interior models by giving a relation between the interior entropy (where the convective zone becomes essentially adiabatic at depth) and the surface parameters $`S(T_{\mathrm{eff}},g)`$. This relation plays a central role in controlling the evolution of fully convective stars. Colors are computed from the synthetic spectra, and are therefore fully consistent with the evolution calculation. This evolution sequence was calculated for objects with solar compositions and masses between 0.01 and 0.3$`M_{}`$, and is very similar to that of Baraffe et al. (1998) since it uses the same input physics (equation of state, atmosphere models, nuclear reaction screening factors, etc.). A limitation of the “NextGen” atmospheres is that they do not include dust opacity, which becomes significant for $`T_{\mathrm{eff}}<2600`$K.

Figure 5 shows the evolution of the absolute magnitudes at $`I`$, $`J`$, $`H`$ and $`K`$ bands, from 1 to 100 MY, based on the models of Saumon & Burrows (unpublished). Each curve shows the evolution for a fixed mass. The two dashed lines highlight the 0.02 and 0.03$`M_{}`$ models. The boxes show the photometric measurements for TWA 5B (Webb et al. 1999, Lowrance et al. 1999), with the height of the box representing the $`\pm 1\sigma `$ photometric error and the width showing the 5–15 MY estimated age of the association. The absolute magnitudes of TWA 5B assume a distance of 55 pc and the $`\pm 9`$pc uncertainty is shown by the error bar in the upper right corner. All four bandpasses indicate that the mass of TWA 5B is between 0.02 and 0.03$`M_{}`$ with an upper limit of $`0.06M_{}`$ if the TWA 5 system lies on the far side of the association and near the upper limit of our age estimate. Given this mass range and the estimated age of TWA 5B, the models indicate that its surface gravity is $`3.8<\mathrm{log}g(\mathrm{cm}/\mathrm{s}^2)<4.0`$.

Stellar and substellar objects with $`M0.012M_{}`$ undergo a phase of nearly constant luminosity which corresponds to the fusion of their primordial deuterium content (D’Antona & Mazzitelli 1985, Saumon et al. 1996). This phase lasts for 2 - 20 MY and contraction — with a consequent steady decrease in luminosity — resumes once the deuterium is exhausted. Figure 5 shows that TWA 5B is almost certainly in the deuterium burning phase of its evolution.

### 4.3 Colors of TWA 5B

The $`IJK`$ colors of TWA 5B are consistent with its dM8.5–dM9 spectral classification based on 0.65–0.75 $`\mu `$m spectra (Webb et al. (1999); Leggett, Allard & Hauschildt (1998)), and thus a temperature of $`T_{\mathrm{eff}}=2600\pm 150`$K (Luhman, Liebert & Rieke (1997); Leggett et al. (1996)). In a $`JH`$ vs. $`HK`$ diagram, TWA 5B falls well outside of the observed sequence of very-low mass stars and brown dwarf candidates in the field (Leggett, Allard & Hauschildt (1998)), while all other members of the association fall along the observed sequence of field stars. This indicates that the $`H`$ magnitude for TWA 5B may be erroneous (by $``$ 1$`\sigma `$) or that its relatively low surface gravity results in a redder $`HK`$ color.

The narrowband infrared colors are shown in Fig. 6 along with the synthetic colors from the “NextGen” spectra. Each curve shows the colors for $`T_{\mathrm{eff}}=2600`$, 2800, 3000 and 3200 K (from left to right) for a fixed gravity. The colors of TWA 5B are shown by the triangle with error bars. For the estimated $`T_{\mathrm{eff}}=2600\pm 150`$K and $`\mathrm{log}g=3.9\pm 0.1`$, there is a reasonable agreement for the F164N$``$F215N color but the models are $`0.4`$ magnitude too blue in F164N$``$F190N. Consequently, TWA 5B is brighter at 1.9 $`\mu `$m than predicted by the models.

The F190N bandpass falls in the middle of a strong H2O absorption band (Fig. 3) whose strength probably is overestimated by the “NextGen” models. Allard et al. (1997) compare a sequence of near-infrared spectra of late M dwarfs with their synthetic spectra. In all cases, the models overestimate the depth of the H2O band, an effect which increases for later spectral types. While an inadequate H2O opacity may be partly responsible for this effect, Tsuji, Ohnaka, & Aoki (1996) have shown that the condensation of dust in atmospheres of low $`T_{\mathrm{eff}}`$ results in a source of continuum opacity which decreases the depth of the water absorption bands. New atmosphere models including dust opacity (Tsuji, Ohnaka, & Aoki; Leggett, Allard & Hauschildt 1998) indicate that its effects on the spectrum (and on broadband colors) become discernible for $`T_{\mathrm{eff}}<2800`$K but remain moderate ($`0.1`$ mag) at the effective temperature of TWA 5B ($`2600`$K). While current models including dust opacity may not fully account for the relatively high F190N flux of TWA 5B, the F164N$``$F190N color of TWA 5B is a strong indication of the presence of dust in its atmosphere.

### 4.4 Astrometry of TWA 5B

At a distance of 55 $`\pm `$ 9 pc, TWA 5B lies at a projected distance from TWA 5A of 108 astronomical units. TWA 5A has a spectral type of M1.5 and a likely mass for the central binary of $``$1.4 M while TWA 5B has an estimated spectral type of M8.5 (Webb et al. (1999)) and a likely mass of $``$25 Mjup. Given this information about the TWA 5 system, the orbital period P of TWA 5B should be P $``$ 1000 years. Thus, the angular motion of TWA 5B, assuming a circular orbit of radius 1$`^{\prime \prime }.`$96 viewed nearly pole-on, would be 0$`^{\prime \prime }.`$013 yr-1 (or 0$`^{\prime \prime }.`$010 yr-1 if TWA 5A is a single star with a mass of 0.7 M and P = 1300 years).

The NICMOS observations of TWA 5 were obtained on 25 April (Lowrance et al. (1999)) and 12 July 1998 (this paper), a difference of $``$0.21 years. In only one-fifth of a year, the orbital motion of TWA 5B would have changed its position relative to TWA 5A by only $``$0$`^{\prime \prime }.`$0027, too small for the positional difference to be measured at these two epochs. Thus, the differences between the positions we measured and the corrected positions from earlier epoch observations reported by Lowrance, Weinberger, & Schneider (1999) are strictly due to the relative accuracies of the different measurements.

Although we have demonstrated that the positional change for TWA 5B is measurement error, not orbital motion, we also have shown that the position reported in this paper is a very accurate “starting” position for TWA 5B. In addition, our results show that it is possible to measure the relative separation of these two objects to an accuracy of only a few thousandths of an arcsec. Thus, the orbital motion of TWA 5B should be measurable to a fairly high degree of accuracy with ground-based observing facilities equipped with adaptive optics or with the refurbished NICMOS camera.

## 5 Summary

To the sensitivity limits of these data, our images reveal no detectable circumstellar disks or infrared reflection nebulae, and no low mass stellar or substellar companions around stars in the five studied TWA systems other than the previously discovered TWA 5B. As for TWA 5B itself, our results suggest that this object has a mass in the range of 0.02–0.03 $`M_{}`$, in good agreement with the work of Lowrance et al. (1999). Finally, while our single epoch observations cannot demonstrate or measure the orbital motion of TWA 5B, they are more than accurate enough to permit the measurement of this motion, in combination with future epoch HST or adaptive optics, ground-based observations, with a baseline of only about a year.

We thank the referee for thoughtful suggestions which improved the clarity of the manuscript, A. Burrows for computing the evolutionary sequences used in this work, and F. Allard, P.H. Hauschildt, I. Baraffe and G. Chabrier for making their synthetic spectra and models available. This research was supported by NSF grant AST93-18970 and NASA grants NAG5-4988 and GO07861.01-96A and is based on observations obtained with the NASA/ESA Hubble Space Telescope at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.

### FIGURE CAPTIONS


text:
# Magnetic properties of frustrated spin ladder

## I Introduction

The frustration of the antiferromagnetic exchange interaction brings about many interesting phenomena in the quantum spin systems, because it generally enhances the quantum fluctuation. It would be valuable to consider the effect of the frustration on the spin ladder, like the materials SrCu2O3 (Ref.), Cu2(C2H12N2)2Cl4 (Refs.) and La6Ca8Cu24O41 (Ref.). They are strongly quantized and have the spin gap. When the next-nearest-neighbor (NNN) exchange interaction appears, the frustration takes place in the system. In the classical limit it is easily shown that the system has two different ordered phases depending on the strength of the NNN exchange and the phase boundary does not change even under external magnetic field. In the quantum system, however, some modifications should exist in the ground state phase diagram, because the spin ladder has no long range order even at $`T=0`$. In this paper, we investigate the frustrated spin ladder by the exact diagonalization of the finite clusters to determine the magnetic phase diagram, even under external field. In addition we consider the possibility of the magnetization plateau , which is predicted by a strong coupling approach.

## II Model and numerical method

The $`S=1/2`$ spin ladder with NNN coupling is described by the Hamiltonian

$``$ $`=`$ $`J_1{\displaystyle \underset{i}{\overset{L}{}}}(𝐒_{1,i}𝐒_{1,i+1}+𝐒_{2,i}𝐒_{2,i+1})`$ (1)
$`+`$ $`J_{}{\displaystyle \underset{i}{\overset{L}{}}}(𝐒_{1,i}𝐒_{2,i})`$ (2)
$`+`$ $`J_2{\displaystyle \underset{i}{\overset{L}{}}}(𝐒_{1,i}𝐒_{2,i+1}+𝐒_{2,i}𝐒_{1,i+1}),`$ (3)
where $`J_1`$, $`J_2`$ and $`J_{}`$ are the coupling constants of the leg, NNN (diagonal) and rung exchange interactions, respectively. We put $`J_{}`$=1 in the following. Using the Lanczos algorithm we numerically solved the ground state of the finite clusters. We also calculated the lowest energy of $``$ for $`_i^L(S_{1,i}^z+S_{2,i}^z)=M`$ , which denotes $`E(M)`$. Using $`E(M)`$, we investigate the magnetic state with $`mM/L`$ under the external field described by $`_Z=H_i^L(S_{1,i}^z+S_{2,i}^z)`$.

## III Two magnetic phases

Consider the nonmagnetic ground state at first. In the classical limit the system has two different ordered phases divided by the first-order phase boundary $`J_2=J_{}/2(=1/2)`$ shown as a dashed line in Fig. 1. $`J_1=1/2`$ is also the boundary because the phase diagram should be symmetric under the exchange of $`J_1`$ and $`J_2`$ (the reflection with respect to the dot-dashed line in Fig. 1).

In the quantum $`S=1/2`$ system we should distinguish the two phases based on the dimer picture; the dimers along the rung and the diagonal, respectively. The former is realized for $`J_2J_{}/2`$, while the latter for $`J_2J_{}/2`$. In the latter phase each two spins coupled by the rung are expected to behave like an effective $`S=1`$ (triplet) object. Thus we call the two phases ‘rung-dimer’ and ‘rung-triplet’, respectively. The phase boundary is easily detected as a level crossing point in the ground state even in small finite clusters. Since the boundary is almost independent of $`L`$, we show only the result of $`L=12`$ as circles in Fig. 1. Our study of the spin correlation function along the rung also supported the above argument and suggested that the boundary is first-order. The results are completely consistent with the recent analysis by the density matrix renormalization group.  (It also indicated the crossover of the phase boundary from first-order to second-order ones for $`J_2<0.287J_1`$, but we don’t consider such a parameter region in this paper.)

Even in the magnetic state under external field the two phases still can be identified by the different canted Néel orders shown in the Figs. 2(a) and (b), respectively, in the classical system. The phase boundary is the same as the nonmagnetic ground state. The quantum system is gapless for $`0<m<1`$ and it might be difficult to distinguish the two phases by the dimer picture. In this case the classical picture is useful because the gapless phase is characterized by the power-law decay of the dominant spin correlation function corresponding to the classical order. Thus the quantum system should also have two phases like the classical limit. The same analysis as the nonmagnetic state indicated the first-order boundary for finite $`m`$. We show the boundaries for $`m`$=1/6, 1/3, 1/2, 2/3 and 5/6 ($`L`$=12) in Fig. 1. They exhibit a small $`m`$ dependence, although it is not so large that a field-induced transition between the two phases is expected to occur in any realistic situations. As the magnetization increases, the boundary tends to approach to the classical limit for most magnetizations. For $`m=1/2`$, however, the boundary exhibits a quite different behavior and it is close to the nonmagnetic one. It implies that the quantum fluctuation is enhanced by the frustration particularly at $`m=1/2`$. Thus we consider the possibility of another spin gap induced by external field, that is observed as a plateau in the magnetization curve at $`m=1/2`$.

## IV Magnetization plateau

We consider the magnetization plateau at $`m=1/2`$. The plateau length $`\mathrm{\Delta }E(M+1)+E(M1)2E(M)`$ is one of useful order parameters to investigate the boundary between the gapless and plateau phases. Since $`\mathrm{\Delta }`$ is the low-lying energy gap, it should obey the relation $`\mathrm{\Delta }1/L`$ in the gapless phase. The scaled plateau $`L\mathrm{\Delta }`$ for several $`L`$ is plotted versus $`J_2`$ with fixed $`J_1`$ to 0.4 in Fig. 3 (a). It suggested that a gapless-gapful transition occurs at $`J_20.2`$. To clarify the feature of the transition, we investigate the central charge $`c`$ of the conformal field theory (CFT) and the critical exponent $`\eta `$. $`\eta `$ is defined by the asymptotic behavior of the spin correlation function $`S_0^+S_r^{}(1)^rr^\eta `$. CFT enables us to estimate $`c`$ and $`\eta `$ from the low-lying energy spectra of finite clusters, using the forms $`E(M)/Lϵ(m)\pi cv_s/6L^2`$ and $`\mathrm{\Delta }\pi v_s\eta /L`$ $`(L\mathrm{})`$, where $`v_s`$ is the sound velocity which is the gradient of the dispersion curve at the origin. After some extrapolation to the infinite length limit, we show the results of $`c`$ and $`\eta `$ for $`J_1=0.4`$ in Fig. 3(b). It justifies that the phase boundary is of the Kosterlitz-Thouless (KT) transition with $`c=1`$ in the gapless phase and $`\eta =1`$ at the critical point. Thus we determine the phase boundary as points with $`\eta =1`$ in the $`J_1`$-$`J_2`$ plane. The result of the KT line is shown as solid symbols in Fig. 4 together with the first-order boundary indicated as open symbols. Fig. 4 is a complete phase diagram at $`m=1/2`$. The plateau phase is surrounded by the KT line and first-order line. The intersection of the two lines is expected to be a tri-critical point. The present analysis suggested that the plateau appears only in the rung-dimer phase. The rung-triplet phase reasonably has no plateau, because it is equivalent to the uniform $`S=1`$ chain.

A necessary condition of the presence of the plateau in general 1D systems was rigorously given by $`Q(Sm)`$. $`Q`$ is the periodicity of the ground state and $`S`$ is the total spin of the unit cell. The present case must hold $`Q=2`$ in the plateau phase at $`m=1/2`$. It suggests that the frustration stabilizes the structure where the singlet and triplet rung bonds are alternating, as is in the case of the zigzag ladder.

Finally we present the magnetization curves for ($`J_1`$,$`J_2`$)= (0.5,0), (0.5,0.3) and (0.5,0.4) in Fig. 5. They were obtained by the size scaling in Ref. applied to the calculated energy spectra of finite systems up to $`L=16`$. The plateau clearly appears at $`m=1/2`$ in the latter two cases.

## V Summary

The antiferromagnetic spin ladder with NNN coupling is investigated by the exact diagonalization of finite clusters. It indicated the existence of the two magnetic phases; the rung-dimer and rung-triplet phases, not only for $`m=0`$ but also in the magnetic state. It is also found that the magnetization plateau possibly appears at $`m=1/2`$ only in the rung-dimer phase.


text:
# On the accuracy and running time of GSAT

## 1 Introduction

The problem of deciding satisfiability of a boolean formula is extensively studied in computer science. It appears prominently, as a prototypical NP-complete problem, in the investigations of computational complexity classes. It is studied by the automated theorem proving community. It is also of substantial interest to the AI community due to its applications in several areas including knowledge representation, diagnosis and planning.

Deciding satisfiability of a boolean formula is an NP-complete problem. Thus, it is unlikely that sound and complete algorithms running in polynomial time exist. However, recent years brought several significant advances. First, fast (although, clearly, still exponential in the worst case) implementations of the celebrated Davis-Putnam procedure \[DP60\] were found. These implementations are able to determine in a matter of seconds the satisfiability of critically constrained CNF formulas with 300 variables and thousands of clauses \[DABC96\]. Second, several fast randomized algorithms were proposed and thoroughly studied \[SLM92, SKC96, SK93, MSG97, Spe96\]. These algorithms randomly generate valuations and then apply some local improvement method in an attempt to reach a satisfying assignment. They are often very fast but they provide no guarantee that, given a satisfiable formula, a satisfying assignment will be found. That is, randomized algorithms, while often fast, are not complete. Still, they were shown to be quite effective and solved several practical large-scale satisfiability problems \[KS92\].

One of the most extensively studied randomized algorithms recently is GSAT \[SLM92\]. GSAT was shown to outperform the Davis-Putnam procedure on randomly generated 3-CNF formulas from the crossover region \[SLM92\]. However, GSAT’s performance on structured formulas (encoding coloring and planning problems) was poorer \[SKC96, SK93, SKC94\]. The basic GSAT algorithm would often become trapped within local minima and never reach a solution. To remedy this, several strategies for escaping from local minima were added to GSAT yielding its variants: GSAT with averaging, GSAT with clause weighting, GSAT with random walk strategy (RWS-GSAT), among others \[SK93, SKC94\]. GSAT with random walk strategy was shown to perform especially well. These studies, while conducted on a wide range of classes of formulas rarely address a critical issue of the likelihood that GSAT will find a satisfying assignment, if one exists, and the running time is studied without a reference to this likelihood. Notable exceptions are \[Spe96\], where RWS-GSAT is compared with a simulated annealing algorithm SASAT, and \[MSG97\], where RSW-GSAT is compared to a tabu search method.

In this paper, we propose a systematic approach for studying the quality of randomized algorithms. To this end, we introduce the concepts of the accuracy and of the running time relative to the accuracy. The accuracy measures how likely it is that a randomized algorithm finds a satisfying assignment, assuming that the input formula is satisfiable. It is clear that the accuracy of GSAT (and any other similar randomized algorithm) grows as a function of time — the longer we let the algorithm run, the better the chance that it will find a satisfying valuation (if one exists). In this paper, we present experimental results that allow us to quantify this intuition and get insights into the rate of growth of the accuracy.

The notion of the running time of a randomized algorithm has not been rigorously studied. First, in most cases, a randomized algorithm has its running time determined by the choice of parameters that specify the number of random guesses, the number of random steps in a local improvement process, etc. Second, in practical applications, randomized algorithms are often used in an interactive way. The algorithm is allowed to run until it finds a solution or the user decides not to wait any more, stops the execution, modifies the parameters of the algorithm or modifies the problem, and tries again. Finally, since randomized algorithms are not complete, they may make errors by not finding satisfying assignments when such assignments exist. Algorithms that are faster may be less accurate and the trade-off must be taken into consideration \[Spe96\].

It all points to the problems that arise when attempting to systematically study the running times of randomized algorithms and extrapolate their asymptotic behavior. In this paper, we define the concept of a running time relative to the accuracy. The relative running time is, intuitively, the time needed by a randomized algorithm to guarantee a postulated accuracy. We show in the paper that the relative running time is a useful performance measure for randomized satisfiability testing algorithms. In particular, we show that the running time of GSAT relative to a prescribed accuracy grows exponentially with the size of the problem.

Related work where the emphasis has been on fine tuning parameter settings \[PW96, GW95\] has shown somewhat different results in regard to the increase in time as the size of the problems grow. The growth shown by \[PW96\] is the retropective variation of maxflips rather than the total number of flips. The number of variables for the 3-CNF randomized instances reported \[GW95\] are $`50,70,100`$. Although our results are also limited by the ability of complete algorithms to determine satisfiable instances, we have results for $`50,100,\mathrm{},400`$ variable instances in the crossover region. The focus in our work is on maintaining accuracy as the size of the problems increase.

Second, we study the dependence of the accuracy and the relative running time on the number of satisfying assignments that the input formula admits. Intuitively, the more satisfying assignments the input formula has, the better the chance that a randomized algorithm finds one of them, and the shorter the time needed to do so. Again, our results quantify these intuitions. We show that the performance of GSAT increases exponentially with the growth in the number of satisfying assignments.

These results have interesting implications for the problem of constructing sets of test cases for experimenting with satisfiability algorithms. It is now commonly accepted that random $`k`$-CNF formulas from the cross-over region are “difficult” from the point of view of deciding their satisfiability. Consequently, they are good candidates for testing satisfiability algorithms. These claims are based on the studies of the performance of the Davis-Putnam procedure. Indeed, on average, it takes the most time to decide satisfiability of CNF formulas randomly generated from the cross-over region. However, the suitability of formulas generated randomly from the cross-over region for the studies of the performance of randomized algorithms is less clear. Our results indicate that the performance of randomized algorithms critically depends on the number of satisfying assignments and much less on the density of the problem. Both under-constrained and over-constrained problems with a small number of satisfying assignments turn out to be hard for randomized algorithms. In the same time, Davis-Putnam procedure, while sensitive to the density, is quite robust with respect to the number of satisfying truth assignments.

On the other hand, there are classes of problems that are “easy” for Davis-Putnam procedure. For instance, Davis-Putnam procedure is very effective in finding 3-colorings of graphs from special classes such as 2-trees (see Section 4 for definitions). Thus, they are not appropriate benchmarks for Davis-Putnam type algorithms. However, a common intuition is that structured problems are “hard” for randomized algorithms \[SKC96, SK93, SKC94\]. In this paper we study this claim for the formulas that encode 3- and 4-coloring problem for 2-trees. We show that GSAT’s running time relative to a given accuracy grows exponentially with the size of a graph. This provides a formal evidence to the “hardness” claim for this class of problems and implies that, while not useful in the studies of complete algorithms such as Davis-Putnam method, they are excellent benchmarks for studying the performance of randomized algorithms.

The main contribution of our paper is not as much a discovery of an unexpected behavior of randomized algorithms for testing satisfiability as it is a proposed methodology for studying them. Our concepts of the accuracy and the relative running time allow us to quantify claims that are often accepted on the basis of intuitive arguments but have not been formally pinpointed.

In the paper, we apply our approach to the algorithm RWS-GSAT from \[SK93, SKC94\]. This algorithm is commonly regarded as one of the best randomized algorithms for satisfiability testing to date. For our experiments we used walksat version 35 downloaded from ftp.research.att.com/dist/ai and run on a SPARC Station 20.

## 2 Accuracy and running time

In this section, we will formally introduce the notion of the accuracy of a randomized algorithm $`𝒜`$. We will then define the concept of the running time relative to accuracy.

Let $``$ be a finite set of satisfiable CNF formulas and let $`𝒫`$ be a probability distribution defined on $``$. Let $`𝒜`$ be a sound algorithm (randomized or not) to test satisfiability. By the accuracy of $`𝒜`$ (relative to the probability space $`(,𝒫)`$), we mean the probability that $`𝒜`$ finds a satisfying assignment for a formula generated from $``$ according to the distribution $`𝒫`$. Clearly, the accuracy of complete algorithms (for all possible spaces of satisfiable formulas) is 1 and, intuitively, the higher the accuracy, the more “complete” is the algorithm for the space $`(,𝒫)`$.

When studying and comparing randomized algorithms that are not complete, accuracy seems to be an important characteristics. It needs to be taken into account — in addition to the running time. Clearly, very fast algorithms that often return no satisfying assignments, even if they exist, are not satisfactory. In fact, most of the work on developing better randomized algorithms can be viewed as aimed at increasing the accuracy of these algorithms. Despite this, the accuracy is rarely explicitly mentioned and studied (see \[Spe96, MSG97\]).

We will propose now an approach through which the running times of randomized satisfiability testing algorithms can be compared. We will restrict our considerations to the class of randomized algorithms designed according to the following general pattern. These algorithms consist of a series of tries. In each try, a truth assignment is randomly generated. This truth assignment is then subject to a series of local improvement steps aimed at, eventually, reaching a satisfying assignment. The maximum number of tries the algorithm will attempt and the length of each try are the parameters of the algorithm. They are usually specified by the user. We will denote by $`MT`$ the maximum number of tries and by $`MF`$ — the maximum number of local improvement steps. Algorithms designed according to this pattern differ, besides possible differences in the values $`MT`$ and $`MF`$, in the specific definition of the local improvement process. A class of algorithms of this structure is quite wide and contains, in particular, the GSAT family of algorithms, as well as algorithms based on the simulated annealing approach.

Let $`A`$ be a randomized algorithm falling into the class described above. Clearly, its average running time on instances from the space $`(,𝒫)`$ of satisfiable formulas depends, to a large degree, on the particular choices for $`MT`$ and $`MF`$. To get an objective measure of the running time, independent of $`MT`$ and $`MF`$, when defining time, we require that a postulated accuracy be met. Formally, let $`a`$, $`0<a1`$, be a real number (a postulated accuracy). Define the running time of $`A`$ relative to accuracy $`a`$, $`t^a`$, to be the minimum time $`t`$ such that for some positive integers $`MT`$ and $`MF`$, the algorithm $`A`$ with the maximum of $`MT`$ tries and with the maximum of $`MF`$ local improvement steps per try satisfies:

1. the average running time on instances from $`(,𝒫)`$ is at most $`t`$, and
2. the accuracy of $`A`$ on $`(,𝒫)`$ is at least $`a`$.

Intuitively, $`t^a`$ is the minimum expected time that guarantees accuracy $`a`$. In Section 3, we describe an experimental approach that can be used to estimate the relative running time.

The concepts of accuracy and accuracy relative to the running time open a number of important (and, undoubtedly, very difficult) theoretical problems. However, in this paper we will focus on an experimental study of accuracy and relative running time for a GSAT-type algorithm. These algorithms follow the following general pattern for the local improvement process. Given a truth assignment, GSAT selects a variable such that after its truth value is flipped (changed to the opposite one) the number of unsatisfied clauses is minimum. Then, the flip is actually made depending on the result of some additional (often again random) procedure.

In our experiments, we used two types of data sets. Data sets of the first type consist of randomly generated 3-CNF formulas \[MSL92\]. Data sets of the second type consist of CNF formulas encoding the $`k`$-colorability problem for randomly generated 2-trees. These two classes of data sets, as well as the results of the experiments, are described in detail in the next two sections.

## 3 Random 3-CNF formulas

Consider a randomly generated 3-CNF formula $`F`$, with $`N`$ variables and the ratio of clauses to variables equal to $`L`$. Intuitively, when $`L`$ increases, the probability that $`F`$ is satisfiable should decrease. It is indeed so \[MSL92\]. What is more surprising, it switches from being close to one to being close to zero very abruptly in a very small range from $`L`$ approximately $`4.25`$ to $`4.3`$. The set of 3-CNF formulas at the cross-over region will be denoted by $`CR(N)`$. Implementations of the Davis-Putnam procedure take, on average, the most time on 3-CNF formulas generated (according to a uniform probability distribution) from the cross-over regions. Thus, these formulas are commonly regarded as good test cases for experimental studies of the performance of satisfiability algorithms \[CA93, Fre96\].

We used seven sets of satisfiable 3-CNF formulas generated from the cross-over regions $`CR(N)`$, $`N=100,150,\mathrm{},400`$. These data sets are denoted by $`DS(N)`$. Each data set $`DS(N)`$ was obtained by generating randomly 3-CNF formulas with $`N`$ variables and $`L=4.30`$ (for $`N=100`$) and $`L=4.25`$ (for $`N150`$) clauses. For each formula, the Davis-Putnam algorithm was then used to decide its satisfiability. The first one thousand satisfiable formulas found in this way were chosen to form the data set.

The random algorithms are often used with much larger values of $`N`$ than we have reported in this paper. The importance of accuracy in this study required that we have only satisfiable formulas (otherwise, the accuracy cannot be reliably estimated). This limited the size of randomly generated 3-CNF formulas used in our study since we had to use a complete satisfiability testing procedure to discard those randomly generated formulas that were not satisfiable. In Section 5, we discuss ways in which hard test cases for randomized algorithms can be generated that are not subject to the size limitation.

For each data set $`DS(N)`$, we determined values for $`MF`$, say $`MF_1,\mathrm{},MF_m`$ and $`MT_1,\mathrm{},`$ $`MT_n`$ for use with RWS-GSAT, big enough to result in the accuracy at least 0.98. For instance, for $`N=100`$, $`MF`$ ranged from $`100`$ to $`1000`$, with the increment of 100, and $`MT`$ ranged from 5 to 50, with the increment of 5. Next, for each combination of $`MF`$ and $`MT`$, we ran RWS-GSAT on all formulas in $`DS(N)`$ and tabulated both the running time and the percentage of problems for which the satisfying assignment was found (this quantity was used as an estimate of the accuracy). These estimates and average running times for the data set $`DS(100)`$ are shown in the tables in Figure 1.

Fixing a required accuracy, say at a level of $`a`$, we then looked for the best time which resulted in this (or higher) accuracy. We used this time as an experimental estimate for $`t^a`$. For instance, there are 12 entries in the accuracy table with accuracy $`0.99`$ or more. The lowest value from the corresponding entries in the running time table is 0.03 sec. and it is used as an estimate for $`t^{0.99}`$.

The relative running times $`t^a`$ for RWS-GSAT run on the data sets $`DS(N)`$, $`N=100,150,\mathrm{},400`$, and for $`a=0.90`$ and $`a=0.95`$, are shown in Figure 2. Both graphs demonstrate exponential growth, with the running time increasing by the factor of 1.5 - 2 for every 50 additional variables in the input problems. Thus, while GSAT outperforms Davis-Putnam procedure for instances generated from the critical regions, if we prescribe the accuracy, it is still exponential and, thus, will quickly reach the limits of its applicability. We did not extend our results beyond formulas with up to 400 variables due to the limitations of the Davis-Putnam procedure, (or any other complete method to test satisfiability). For problems of this size, GSAT is still extremely effective (takes only about 2.5 seconds). Data sets used in Section 5 do not have this limitation (we know all formulas in these sets are satisfiable and there is no need to refer to complete satisfiability testing programs). The results presented there also illustrate the exponential growth of the relative running time and are consistent with those discussed here.

## 4 Number of satisfying assignments

It seems intuitive that accuracy and running time would be dependent on the number of possible satisfying assignments. Studies using randomly generated 3-CNF formulas \[CFG<sup>+</sup>96\] and 3-CNF formulas generated randomly with parameters allowing the user to control the number of satisfiable solutions for each instance \[CI95\] show this correlation.

In the same way as for the data sets $`DS(N)`$, we constructed data sets $`DS(100,p_{k1},p_k)`$, where $`p_0=1`$, and $`p_k=2^{k3}100`$, $`k=2,\mathrm{},11`$. Each data set $`DS(100,p_{k1},p_k)`$ consists of 100 satisfiable 3-CNF formulas generated from the cross-over region $`CR(100)`$ and having more than $`p_{k1}`$ and no more than $`p_k`$ satisfying assignments. Each data set was formed by randomly generating 3-CNF formulas from the cross-over region $`CR(100)`$ and by selecting the first 100 formulas with the number of satisfying assignments falling in the prescribed range (again, we used the Davis-Putnam procedure, here).

For each data set we ran the RWS-GSAT algorithm with $`MF=500`$ and $`MT=50`$ thus, allowing the same upper limits for the number of random steps for all data sets (these values resulted in the accuracy of .99 in our experiments with the data set $`DS(100)`$ discussed earlier). Figure 3 summarizes our findings. It shows that there is a strong relationship between accuracy and the number of possible satisfying assignments. Generally, instances with small number of solutions are much harder for RWS-GSAT than those with large numbers of solutions. Moreover, this observation is not affected by how constrained the input formulas are. We observed the same general behavior when we repeated the experiment for data sets of 3-CNF formulas generated from the under-constrained region (100 variables, 410 clauses) and over-constrained region (100 variables, 450 clauses), with under-constrained instances with few solutions being the hardest.

These results indicate that, when generating data sets for experimental studies of randomized algorithms, it is more important to ensure that they have few solutions rather than that they come from the critically constrained region.

## 5 CNF formulas encoding $`k`$-colorability

To expand the scope of applicability of our results and argue their robustness, we also used in our study data sets consisting of CNF formulas encoding the $`k`$-colorability problem for graphs. While easy for Davis-Putnam procedure (which resolves their satisfiability in polynomial time), formulas of this type are believed to be “hard” for randomized algorithms and were used in the past in the experimental studies of their performance. In particular, it was reported in \[SK93\] that RWS-GSAT does not perform well on such inputs (see also \[JAMS91\]).

Given a graph $`G`$ with the vertex set $`V=\{v_1,\mathrm{},v_n\}`$ and the edge set $`E=\{e_1,\mathrm{},e_m\}`$, we construct the CNF formula $`COL(G,k)`$ as follows. First, we introduce new propositional variables $`col(v,i)`$, $`vV`$ and $`i=1,\mathrm{},k`$. The variable $`col(v,i)`$ expresses the fact that the vertex $`v`$ is colored with the color $`i`$. Now, we define $`COL(G,k)`$ to consist of the following clauses:

1. $`\neg col(x,i)\neg col(y,i)`$, for every edge $`\{x,y\}`$ from $`G`$,
2. $`col(x,1)\mathrm{}col(x,k)`$, for every vertex $`x`$ of $`G`$,
3. $`\neg col(x,i)\neg col(x,j)`$, for every vertex $`x`$ of $`G`$ and for every $`i,j`$, $`1i<jk`$.

It is easy to see that there is a one-to-one correspondence between $`k`$-colorings of $`G`$ and satisfying assignments for $`COL(k,G)`$. To generate formulas for experimenting with RWS-GSAT (and other satisfiability testing procedures) it is, then, enough to generate graphs $`G`$ and produce formulas $`COL(G,k)`$.

In our experiments, we used formulas that encode $`3`$-colorings for graphs known as $`2`$-trees. The class of 2-trees is defined inductively as follows:

1. A complete graph on three vertices (a “triangle”) is a 2-tree
2. If $`T`$ is a 2-tree than a graph obtained by selecting an edge $`\{x,y\}`$ in $`T`$, adding to $`T`$ a new vertex $`z`$ and joining $`z`$ to $`x`$ and $`y`$ is also a 2-tree.

A 2-tree with 6 vertices is shown in Fig. 4. The vertices of the original triangle are labeled 1, 2 and 3. The remaining vertices are labeled according to the order they were added.

The concept of 2-trees can be generalized to $`k`$-trees, for an arbitrary $`k2`$. Graphs in these classes are important. They have bounded tree-width and, consequently, many NP-complete problems can be solved for them in polynomial time \[AP89\].

We can generate 2-trees randomly by simulating the definition given above and by selecting an edge for “expansion” randomly in the current 2-tree $`T`$. We generated in this way families $`G(p)`$, for $`p=50,60,\mathrm{},150`$, each consisting of one hundred randomly generated 2-trees with $`p`$ vertices. Then, we created sets of CNF formulas $`C(p,3)=\{COL(T,3):TG(p)\}`$, for $`p=50,60,\mathrm{},150`$. Each formula in a set $`C(p,3)`$ has exactly 6 satisfying assignments (since each 2-tree has exactly 6 different 3-colorings). Thus, they are appropriate for testing the accuracy of RWS-GSAT.

Using CNF formulas of this type has an important benefit. Data sets can be prepared without the need to use complete (but very inefficient for large inputs) satisfiability testing procedures. By appropriately choosing the underlying graphs, we can guarantee the satisfiability of the resulting formulas and, often, we also have some control over the number of solutions (for instance, in the case of 3-colorability of 2-trees there are exactly 6 solutions).

We used the same methodology as the one described in the previous section to tabulate the accuracy and the running time of RSW-GSAT for a large range of choices for the parameters $`MF`$ and $`MT`$. Based on these tables, as before, we computed estimates for the times $`t^a`$ for $`a=0.95`$, for each of the data sets. The results that present the running time $`t^a`$ as a function of the number of vertices in a graph (which is of the same order as the number of variables in the corresponding CNF formula) are gathered in Figure 5. They show that RWS-GSAT’s performance deteriorates exponentially (time grows by the factor of $`34`$ for every 50 additional vertices).

An important question is: how to approach constraint satisfaction problems if they seem to be beyond the scope of applicability of randomized algorithms? A common approach is to relax some constraints. It often works because the resulting constraint sets (theories) are “easier” to satisfy (admit more satisfying assignments). We have already discussed the issue of the number of solutions in the previous section. Now, we will illustrate the effect of increasing the number of solutions (relaxing the constraints) in the case of the colorability problem. To this end, we will consider formulas from the spaces $`C(p,4)`$, representing 4-colorability of 2-trees. These formulas have exponentially many satisfying truth assignments (a 2-tree with $`p`$ vertices has exactly $`3\times 2^p`$ 4-colorings). For these formulas we also tabulated the times $`t^a`$, for $`a=0.95`$, as a function of the number of vertices in the graph. The results are shown in Figure 6.

Thus, despite the fact the size of a formula from $`C(p,4)`$ is larger than the size of a formula from $`C(p,3)`$ by the factor of $`1.6`$, RWS-GSAT’s running times are much lower. In particular, within .5 seconds RWS-GSAT can find a 4-coloring of randomly generated 2-trees with 500 vertices. As demonstrated by Figure 5, RWS-GSAT would require thousands of seconds for 2-trees of this size to guarantee the same accuracy when finding 3-colorings. Thus, even a rather modest relaxation of constraints can increase the number of satisfying assignments substantially enough to lead to noticeable speed-ups. On the other hand, even though “easier”, the theories encoding the 4-colorability problem for 2-trees still are hard to solve by GSAT as the rate of growth of the relative running time is exponential (Fig. 6).

The results of this section further confirm and provide quantitative insights into our earlier claims about the exponential behavior of the relative running time for GSAT and on the dependence of the relative running time on the number of solutions. However, they also point out that by selecting a class of graphs (we selected the class of 2-trees here but there are, clearly, many other possibilities) and a graph problem (we focused on colorability but there are many other problems such as hamiltonicity, existence of vertex covers, cliques, etc.) then encoding these problems for graphs from the selected class yields a family of formulas that can be used in testing satisfiability algorithms. The main benefit of the approach is that by selecting a suitable class of graphs, we can guarantee satisfiability of the resulting formulas and can control the number of solutions, thus eliminating the need to resort to complete satisfiability procedures when preparing the test cases. We intend to further pursue this direction.

## 6 Conclusions

In the paper we formally stated the definitions of the accuracy of a randomized algorithm and of its running time relative to a prescribed accuracy. We showed that these notions enable objective studies and comparisons of the performance and quality of randomized algorithms. We applied our approach to study the RSW-GSAT algorithm. We showed that, given a prescribed accuracy, the running time of RWS-GSAT was exponential in the number of variables for several classes of randomly generated CNF formulas. We also showed that the accuracy (and, consequently, the running time relative to the accuracy) strongly depended on the number of satisfying assignments: the bigger this number, the easier was the problem for RWS-GSAT. This observation is independent of the “density” of the input formula. The results suggest that satisfiable CNF formulas with few satisfying assignments are hard for RWS-GSAT and should be used for comparisons and benchmarking. One such class of formulas, CNF encodings of the 3-colorability problem for 2-trees was described in the paper and used in our study of RWS-GSAT.

Exponential behavior of RWS-GSAT points to the limitations of randomized algorithms. However, our results indicating that input formulas with more solutions are “easier” for RWS-GSAT to deal with, explain RWS-GSAT’s success in solving some large practical problems. They can be made “easy” for RWS-GSAT by relaxing some of the constraints.


text:
# Introduction

## Introduction

The Standard Model (SM) is very successful in confronting the data coming from the highest energy accelerators. Still, there are theoretical reasons to expect that it is not complete, and one of the first questions in the quest for new physics is what is the relevant scale, where new phenomena can give experimental signatures. Recently, a radical proposal  has been put forward for the solution of the hierarchy problem, which brings close the electroweak scale $`\mathrm{m}_{\mathrm{EW}}1\mathrm{TeV}`$ and the Planck scale $`\mathrm{M}_{\mathrm{Pl}}=\frac{1}{\sqrt{\mathrm{G}_\mathrm{N}}}10^{15}\mathrm{TeV}`$. In this framework the effective four-dimensional $`\mathrm{M}_{\mathrm{Pl}}`$ is connected to a new $`\mathrm{M}_{\mathrm{Pl}(4+\mathrm{n})}`$ scale in a (4+n) dimensional theory:

$$\mathrm{M}_{\mathrm{Pl}}^2\mathrm{M}_{\mathrm{Pl}(4+\mathrm{n})}^{2+\mathrm{n}}\mathrm{R}^\mathrm{n}$$
(1)
where there are n extra compact spatial dimensions of radius $`\mathrm{R}`$. This can explain the observed weakness of gravity at large distances. At the same time, quantum gravity becomes strong at a scale M of the order of 1 TeV and could have observable signatures at present and future colliders.

The first experimental searches for large extra dimensions have concentrated on the effects of real and virtual graviton emission111For searches in Bhabha scattering see e.g. .. In a string theory of quantum gravity  there are additional modifications of Standard Model amplitudes and new phenomenological consequences. Effective contact interactions caused by massive string mode oscillations might compete with or even become stronger than those due to virtual exchange of Kaluza-Klein excitations of gravitons, and thus provide the first signature of low scale gravity or a lower bound on the string scale.

Bhabha scattering above the Z resonance offers a reach hunting field for new phenomena . It can be used to search for manifestations of contact interactions and as a very sensitive probe of the point-like structure of the electron.

This paper is organized as follows. In sections 2 and 3 the experimental data and the analysis technique are presented. In the following section, we describe the search for effects of TeV strings in Bhabha scattering. In sections 5 and 6 we use the data to obtain limits on the scale of different contact interaction models, and on the size of electrons respectively. We conclude with a discussion of the results.

## Experimental Data

Data on fermion-pair production at 183 or 189 GeV centre-of-mass energies from the LEP2 collider has become available recently. In the following we will concentrate on the measurements of Bhabha scattering at these two highest energy points, where large data samples have been accumulated during the very successful LEP runs in 1997 and 1998.

The ALEPH , L3  and OPAL  collaborations have presented results for the differential cross section of Bhabha scattering. In the case of L3 and OPAL the results are for both energy points and the scattering angle $`\theta `$ is the angle between the incoming and the outgoing electrons in the laboratory frame. In the ALEPH case the measurements are at 183 GeV and the scattering angle is defined in the outgoing $`\mathrm{e}^+\mathrm{e}^{}`$ rest frame. The acceptance is given by the angular range $`|\mathrm{cos}\theta |<0.9`$ for the ALEPH and OPAL measurements and by $`44^{}<\theta <136^{}`$ for the L3 measurement.

The experiments use different strategies to isolate the high energy sample, where the interactions take place at energies close to the full available centre-of-mass energy. This sample is the main search field for new physics. L3 and OPAL apply an acollinearity cut of $`25^{}`$ and $`10^{}`$ respectively. ALEPH defines the effective energy, $`\sqrt{\mathrm{s}^{}}`$, as the invariant mass of the outgoing fermion pair. It is determined from the angles of the outgoing fermions. For details of the selection procedures, the statistical and systematic errors we refer the reader to the publications of the LEP experiments.

## Analysis Method

The Standard Model predictions for the differential cross sections of Bhabha scattering at 183 and 189 GeV are computed with the Monte Carlo generator BHWIDE . We assign a theory uncertainty of 1.5 % on the absolute scale of the predictions. In all cases the individual experimental cuts of the selection procedures and the isolation of the high energy samples are taken into account. The results are cross-checked with the semi-analytic program TOPAZ0 .

The effects of new phenomena are computed as a function of a generic parameter $`\epsilon `$, defined for each individual case in the corresponding section. Initial-state radiation (ISR) changes the effective centre-of-mass energy in a large fraction of the observed events. We take these effects into account by computing the first order exponentiated differential cross section following . Other QED and electroweak corrections give smaller effects and are neglected.

In total we have 47 data points: 28 from the 3 differential spectra at 183 GeV and 19 from the L3 and OPAL spectra at 189 GeV. A fitting procedure similar to the one in  is applied.

A negative log-likelihood function is constructed by combining all data points at the two centre-of-mass energies:

$$\mathrm{log}=\underset{\mathrm{r}=1}{\overset{\mathrm{n}}{}}\left(\frac{(\mathrm{Prediction}(\mathrm{SM},\epsilon )\mathrm{Measurement})^2}{2\mathrm{\Delta }_{\mathrm{Measurement}}^2}\right)_\mathrm{r}$$
(2)
$`\mathrm{\Delta }_{\mathrm{Measurement}}`$ $`=`$ $`\mathrm{error}(\mathrm{Prediction}(\mathrm{SM},\epsilon )\mathrm{Measurement})`$ (3)
where $`Prediction(SM,\epsilon )`$ is the SM expectation for a given measurement (a point in the differential spectra) combined with the additional effect of new phenomena as a function of the mass scale or electron size, and $`Measurement`$ is the corresponding measured quantity. The index $`\mathrm{r}`$ runs over all data points. The error on a deviation consists of three parts, which are combined in quadrature: a statistical error and a systematic error (as given by the experiments) and the theoretical error assigned above. The systematic errors account for small correlations between data points.

## TeV Strings in Bhabha Scattering

In  the authors develop a model to study the effects of string Regge excitations on physical cross sections by a simple embedding of the Quantum Electrodynamics of electrons and photons into string theory. They use only one gauge group and only vector-like couplings, in order to avoid complications but grasp the general phenomenological picture. The results are model-dependent.

The effects of TeV scale strings on Bhabha scattering are computed from the leading-order scattering amplitudes. All amplitudes are multiplied by a common form-factor

$$𝒮(\mathrm{s},\mathrm{t})=\frac{\mathrm{\Gamma }(1\frac{\mathrm{s}}{\mathrm{M}_\mathrm{S}^2})\mathrm{\Gamma }(1\frac{\mathrm{t}}{\mathrm{M}_\mathrm{S}^2})}{\mathrm{\Gamma }(1\frac{\mathrm{s}}{\mathrm{M}_\mathrm{S}^2}\frac{\mathrm{t}}{\mathrm{M}_\mathrm{S}^2})}.$$
(4)
In the case where the string scale $`\mathrm{M}_\mathrm{S}`$ is close to or smaller than the centre-of-mass energy, the Gamma-functions in this form-factor produce a very reach and complicated resonance structure. On the other hand, in the limit where the Mandelstam variables s and t are much smaller than $`\mathrm{M}_\mathrm{S}`$, we have

$$𝒮(\mathrm{s},\mathrm{t})=(1\frac{\pi ^2}{6}\frac{\mathrm{st}}{\mathrm{M}_\mathrm{S}^4}+\mathrm{}).$$
(5)
So in this model the leading corrections are proportional to $`\mathrm{M}_\mathrm{S}^4`$, corresponding to an operator of dimension 8.

To compare the string predictions to the data on Bhabha scattering above the Z resonance one has to handle also the contributions due to Z exchange and the interference with photon exchange amplitudes. The Z is not part of the string QED model developed in , but as all QED Bhabha scattering amplitudes are multiplied by the common factor $`𝒮(\mathrm{s},\mathrm{t})`$, the authors suggest to compare the differential cross section to the simple formula

$$\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{cos}\theta }=(\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{cos}\theta })_{\mathrm{SM}}|𝒮(\mathrm{s},\mathrm{t})|^2.$$
(6)
The data from the LEP collaborations at 183 and 189 GeV show no statistically significant deviations from the SM predictions due to string effects. In their absence, we use the log-likelihood method, which after proper normalization gives the confidence level for any value of the scale $`\mathrm{M}_\mathrm{S}`$ in the physically allowed region. The exact definition can be found in . The one-sided lower limit on the scale $`\mathrm{M}_\mathrm{S}`$ at 95% confidence level is:

$$\mathrm{M}_\mathrm{S}=0.631\mathrm{TeV}.$$
(7)
Examples of the data analysis at 189 GeV are shown in Figure 1 and Figure 2, where the SM predictions and the expectations from several manifestations of new phenomena are compared to the measurements of the L3 and OPAL collaborations, respectively. In these figures we plot the combined statistical and systematic errors; the theory uncertainty is not shown. In the area of the forward peak the theory uncertainty in the SM prediction starts to limit the precision of our study.

## Contact Interactions

The standard framework, used in searches for deviations from the SM predictions, is the most general combination of helicity conserving dimension-6 operators . In this scheme, new interactions beyond the Standard Model are characterised by a coupling strength, $`g`$, and by an energy scale, $`\mathrm{\Lambda }`$, which can be viewed as the scale of compositeness. At energies much lower than $`\mathrm{\Lambda }`$, we have an effective Lagrangian leading to four-fermion contact interactions.

The differential cross section for fermion-pair production in $`\mathrm{e}^+\mathrm{e}^{}`$ collisions can be decomposed in the usual way as:

$$\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }}=\mathrm{SM}(\mathrm{s},\mathrm{t})+\epsilon \mathrm{C}_{\mathrm{Int}}(\mathrm{s},\mathrm{t})+\epsilon ^2\mathrm{C}_{\mathrm{CI}}(\mathrm{s},\mathrm{t})$$
(8)
where $`\mathrm{SM}(\mathrm{s},\mathrm{t})`$ is the Standard Model contribution, $`\mathrm{C}_{\mathrm{CI}}(\mathrm{s},\mathrm{t})`$ comes from the contact interaction amplitude and $`\mathrm{C}_{\mathrm{Int}}(\mathrm{s},\mathrm{t})`$ is the interference between the SM and the contact interaction terms. The exact form of these functions is given in . By convention $`\frac{\mathrm{g}^2}{4\pi }=1`$ and $`|\eta _{\mathrm{ij}}|1`$, where $`(\mathrm{i},\mathrm{j}=\mathrm{L},\mathrm{R})`$ labels the helicity of the incoming and outgoing fermions. We define

$$\epsilon =\frac{\mathrm{g}^2}{4\pi }\frac{\mathrm{sign}(\eta )}{\mathrm{\Lambda }^2}$$
(9)
where the sign of $`\eta `$ enables to study both the cases of positive and negative interference.

As discussed in the previous section, the data from the LEP collaborations at 183 and 189 GeV show no statistically significant deviations from the SM predictions. In their absence, using the same technique we derive one-sided lower limits on the scale $`\mathrm{\Lambda }`$ of contact interactions at 95% confidence level. They are summarized in Table 1 and Figure 3. The results presented here improve on the limits obtained by individual LEP experiments .

## Electron Size

In the Standard Model leptons, quarks and gauge bosons are considered as point-like particles. A possible substructure or new interactions at as yet unexplored very high energies could manifest themselves as finite radii and anomalous magnetic dipole moments of these particles.

The high precision measurements of the magnetic dipole moment $`(\mathrm{g}2)_\mathrm{e}`$ of the electron can be used to put stringent limits on the electron radius $`\mathrm{r}_\mathrm{e}`$ . If non-standard contributions to $`(\mathrm{g}2)_\mathrm{e}`$ scale linearly with the electron mass, the bound is $`\mathrm{r}_\mathrm{e}210^{23}\mathrm{m}`$. On the other hand, if they scale quadratically with the electron mass, which is a natural consequence of chiral symmetry , the bound is reduced to $`\mathrm{r}_\mathrm{e}310^{18}\mathrm{m}`$. In  the authors perform an analysis of the high precision data on the Z resonance, noting that while the assumption of elementary photons is quite natural, the same is less obvious for the very massive Z bosons. In the pure electron case the limit is not competitive with the $`(\mathrm{g}2)_\mathrm{e}`$ results.

Here we perform a new analysis based on the LEP2 data on Bhabha scattering, where again the photon exchange gives the dominating amplitudes both in the t- and s-channels, and good sensitivity to electron substructure can be expected. The differential cross section for fermion-pair production in $`\mathrm{e}^+\mathrm{e}^{}`$ collisions far above the Z is modified as:

$$\frac{\mathrm{d}\sigma }{\mathrm{dQ}^2}=(\frac{\mathrm{d}\sigma }{\mathrm{dQ}^2})_{\mathrm{SM}}\mathrm{F}_\mathrm{e}^2(\mathrm{Q}^2)\mathrm{F}_\mathrm{f}^2(\mathrm{Q}^2)$$
(10)
where $`\mathrm{F}_\mathrm{e}`$ and $`\mathrm{F}_\mathrm{f}`$ are the form-factors of the initial (final) state fermions. They are parametrized in the standard way as :

$$\mathrm{F}(\mathrm{Q}^2)=1+\frac{1}{6}\mathrm{Q}^2\mathrm{r}^2$$
(11)
where $`\mathrm{Q}^2`$ is the Mandelstam variable s or t for s- or t-channel exchange, and $`\mathrm{r}^2`$ is the mean-square radius of the fermions. This formalism is a convenient way to estimate the electron size in the case where the product $`\mathrm{Q}^2\mathrm{r}^2`$ is small.

From the data of the LEP collaborations at 183 and 189 GeV we extract the following upper limit on the electron radius at 95% confidence level:

$$\mathrm{r}_\mathrm{e}<2.810^{19}\mathrm{m}.$$
(12)
This limit is one order of magnitude lower than the limit derived from $`(\mathrm{g}2)_\mathrm{e}`$ measurements in the case where the deviations from the SM of the magnetic dipole moment of the electron depend quadratically on its mass.

High energy analyses have been performed in interactions involving electrons and quarks, assuming a single form-factor for all fermions. The H1 collaboration at HERA uses deep inelastic scattering and obtains a limit of $`\mathrm{r}<2610^{19}\mathrm{m}`$ at 95 % confidence level . The CDF collaboration at the TEVATRON studies the Drell-Yan process to put a limit of $`\mathrm{r}<5.610^{19}\mathrm{m}`$ at 95 % confidence level .

## Discussion

The search for TeV strings motivates a fresh look at Bhabha scattering. In the model analyzed here the string realization of quantum gravity is manifested as a form-factor which modifies the differential cross section. The lower limit obtained in our analysis of LEP2 data is $`\mathrm{M}_\mathrm{S}=0.631\mathrm{TeV}`$. In  from the study of virtual graviton exchange in gravity models with large extra dimensions we obtained a lower limit on their scale of $`\mathrm{\Lambda }_\mathrm{T}=1.412\mathrm{TeV}`$ for positive inteference ($`\lambda =+1`$222 This value of $`\mathrm{\Lambda }_\mathrm{T}`$ corresponds, depending on the convention, also to a gravity scale $`\mathrm{M}_\mathrm{s}=1.261\mathrm{TeV}`$. The gravity scale with subscript small s should not be confused with the string scale $`\mathrm{M}_\mathrm{S}`$, studied here.. As noted in , the gravity scale is between $`1.6÷3.0\mathrm{M}_\mathrm{S}`$, depending on the coupling strength. The results on the gravity scale from  and on the string scale from this analysis agree well with each other.

It is interesting to note that our study of the electron size also leads to form-factors modifying the differential cross section, but with opposite sign. The limit derived here, $`\mathrm{r}_\mathrm{e}<2.810^{19}\mathrm{m}`$, becomes $`\mathrm{M}_\mathrm{r}>0.705\mathrm{TeV}`$, if translated to a mass scale. This is a reflection of the similar magnitude of the effects at LEP2 energies in both cases, even if the physics mechanisms involved are different.

In the framework of contact interactions very stringent bounds exceeding 10 TeV are obtained. When interpreting the physical meaning of these limits, we should remember that a strong coupling $`\frac{\mathrm{g}^2}{4\pi }=1`$ for the novel interactions is postulated by convention. If we assume a coupling of electromagnetic strength, the limits can be translated:

$$\mathrm{\Lambda }^{}=\sqrt{\alpha _{\mathrm{QED}}}\mathrm{\Lambda }=0.085\mathrm{\Lambda }$$
(13)
where we have used the value of the fine structure constant and ignored the small effect of a running $`\alpha _{\mathrm{QED}}`$. For instance the VV model with positive interference gives effects similar to the ones resulting from a finite electron size, as shown in Figure 1 and Figure 2. The limit for the VV model translates as follows:

$$\mathrm{\Lambda }_+=13.0\mathrm{TeV}\mathrm{\Lambda }^{}=1.1\mathrm{TeV}\mathrm{r}=1.810^{19}\mathrm{m}.$$
(14)
This results is comparable with the upper limit for electron substructure, derived using form-factors.

The measurements of Bhabha scattering above the Z resonance confirm the predictions of the Standard Model and reach already a similar level of precision as the best theoretical tools available. In order to fully exploit the physics potential of the large data samples collected during the LEP running in 1999 and expected in 2000, improved theory predictions are very desirable. Bhabha scattering is a probe, sensitive enough to provide a first window to new physics at the TeV scale.

## Acknowledgements

The author is grateful to A. Böhm, M. Peskin and I. Antoniadis for valuable discussions.