text:
# KMC-MD Investigations of Hyperthermal Copper Deposition on Cu(111)

## I Introduction

Recent technological advances have resulted in wide-ranging implementation of devices which utilize hyperthermal energy particles for thin-film and nano-scale device growthMinv . While these advances have propelled production efforts, a detailed understanding of the relevant physical mechanisms has not been fully developed. It has become increasingly apparent that even hyperthermal energetic particles can stimulate a variety of thermal and non-thermal processes, ranging from strain relaxationKarr  to sputter erosionMich2 , enhanced nuclei densitiesEsch , and improved composite layer adhesionFlad1 .

Efforts to model crystal surfaces during deposition have provided detailed information about the nature of surface kinetics. For example, molecular dynamics has allowed accurate predictions of many energy barriers for surface diffusion processesStolze . Many of these predictions have been carefully addressed experimentallyMorg1 ; Rosen3 ; Geis1 , but most efforts use deposition techniques with thermally generated constituents arriving at the substrate with less than a tenth of an eV.

Many of these studies have identified the “Ehrlich-Schwöebel” ehrl1 ; schw1  barrier for interlayer diffusion as the source of three dimensional growth modes. The resulting “uphill” current produced by this interlayer diffusion barrier can be reduced by increasing temperature to surmount this barrierRosen1 ; Wulf1  or sufficiently decreasing step-step separationgeis2 . However, in hetero-structures and nano-structures, elevated temperatures result in inter-diffusion, chemical reaction, and thermodynamic relaxation, making nano-scale patterns difficult to retain.

Physical vapor deposition (PVD) techniques utilizing hyperthermal energy constituents produce smoother epitaxial films in many systems with finite Ehrlich-Schwöebel barriersEsch2 ; Got1 . It has been proposed that a “peening” effect known for building-up stress at medium energiesSun ; Wind1 ; Wind2  may act as a relaxation mechanism at low and hyperthermal energiesKoster .

Progress toward understanding hyperthermal energy collisions has been hindered by a lack of models which both accurately describe the collision process and the kinetic processes at realistic deposition rates ($`1`$ ML/s). The KMC-MD method allows the complexity of the atomic collision to be modeled uniquely with molecular dynamics for each atom without prior bias. Between deposition events, the kinetic Monte Carlo evolves the system using well understood kinetics, until the next deposition event. Previous KMC-MD studies of platinum and silver were able to provided new insights into the role of the hyperthermal atom collision during growthJJ .

This paper presents results for energetic collisions on the Cu(111) surface: first, isolated molecular dynamics studies of atom impacts, and, second, KMC-MD simulations of entire films grown with hyperthermal energy atoms. A brief description of the simulation, improvements, and the energy barriers used for our simulations is presented first. The yields for various atomistic mechanisms resulting from molecular dynamics simulations of isolated impacts is presented next. The results of the full KMC-MD growth model at 50 K are then discussed. Finally, simulations of sub-monolayer films at room temperature are presented.

## II Simulation Details

Accurate modeling of crystal growth with hyperthermal energy atoms requires accurately modeling two classes of events active at times scales separated by about six orders of magnitude. The hyperthermal atom impact and subsequent thermalization process is complete in about four picoseconds. On the other hand, the surface kinetics are active on the microsecond time scale. The technical challenge of accurate modeling sub-picosecond events for microseconds has been resolved by coupling two techniques, one appropriate for each time scale. Molecular dynamics is used to simulate the hyperthermal atom collisions, but is not feasible for modeling surface diffusion at realistic deposition rates. For surface diffusion we use kinetic Monte Carlo, which passes atomic configurations at selected impact sites to the molecular dynamics, and accepts the new configurations following an impact. The full details of these simulations have been detailed elsewhere JJ .

We use a bowl-shaped molecular dynamics cluster with three classes of atoms: fully dynamic atoms nearest to the site of the impact, surrounded by three layers of dynamic Langevin atoms, and finally, an outer shell with four layers of static atoms. The Langevin coefficients and cell size are tuned to dissipate energy which can reflect from the boundary and lead to unrealistic rates for energetic reflection and adatom/vacancy formation.

We have found it useful to model isolated hyperthermal atom collisions with a few selected atomic configurations to develop a general picture of which atomic mechanisms are important at various energies. Many atomic configurations on the surface can be classified according to the distance to a step edge, so we have selected several positions near a (111) step edge. We model 100 collisions in each of five atomic cells above a step edge, at the step edge, and five cells below the step edge for each energy of interest. For each collision, impact parameters are randomly selected, and the cluster evolved until the Langevin atoms have thermalized the system. The configuration is then frozen into a final state for analysis and saved. Once satisfactory statistics are developed at one position, the impact site is moved one atomic cell, and the process is repeated. We have found that the statistics for the fifth cell above or below is representative of all cells further from the step.

The role of the molecular dynamics during KMC-MD is the same, except that the configuration of atoms in the molecular dynamics cell is determined by the local environment of the impact site. The KMC uses a hexagonal lattice with in-plane periodic boundary conditions and a compliment of 23 pre-defined thermal moves. Interstitial and HCP lattice positions are not allowed. The activation energies for these moves were calculated using ARTworkartwork  and are listed in Table 1. (For the 273 K simulations the barriers in Table 2 are included as well.) Included in the KMC rate table is a flux weighted choice for adding new atoms. When the algorithm chooses to introduce a new atom, an impact site is randomly selected. The local configuration of atoms is then copied into the molecular dynamics cluster, which simulates the collision. Once the final configuration is determined, it is returned to the KMC, and thermal evolution continues.

### II.1 Modifications to the KMC-MD algorithm

One of the trickier parts of the KMC-MD method involves moving atoms from the continuous MD space to the discrete lattice of the KMC. During collisions involving many atoms on the surface, clusters of atoms occasionally freeze into HCP rather than FCC lattice positions. (HCP and FCC lattice positions are energetically equivalent on the (111) face using the EMT potential for copper.) In the previous studies of platinum and silverJJ , these clusters rarely exceeded five atoms. As the algorithm encountered atoms in HCP sites it would place them on the nearest available FCC site. During simulations of copper, these clusters are sometimes as large as eight atoms, and this previous technique did not always preserve the shape of the cluster. In some cases, atoms near the middle of a cluster could not be placed at all, since all nearby sites would already be filled.

Our modified algorithm creates a list of atoms in HCP sites, then sorts the list from highest coordination to lowest. As each atom in the list is selected for placement, the three fcc sites surrounding that HCP site are checked for occupancy. The wayward atom is then placed into the unoccupied fcc site with the highest coordination (a random selection occurs if multiple sites have the highest coordination.) Since atoms near the center of a cluster get placed first, all atoms have an available FCC position. This change preserves the cluster and has successfully placed all the atoms in supported FCC sites.

## III Results for MD Collisions

Simulating individual atomic collisions in a pre-selected environment can provide a general insight into the yields for atomic mechanisms at different energies. Once an impact has been simulated, the final atomic configuration is classified according to the change in the population of the atomic layers. If the impact site is above the step edge, and the atom incident atom is incorporated into the step, the event is considered an insertion. We do not distinguish between an atom that is actually inserted, and one which just bounces over the step edge. A decrease in the population of a layer requires the formation of a vacancy. The formation of vacancies usually provides additional adatoms (adatom-vacancy pairs) that can contribute to surface relaxation through enhanced lateral diffusion. If the total number of atoms in the cluster decreases, this is considered a sputtering event. (Spontaneous thermodynamic re-evaporation is negligible.)

The yields averaged from many simulated collisions at several different atomic configurations are presented in Figure 1. The hyperthermal atomic mechanisms observed in platinum and silver are present in copper, but the specific energies of activation varying somewhat. The insertion mechanism is active at the first position above the step edge as low as 3 eV. As the energy increases, atoms are inserted deeper into the island, which increases the yield. By 9 eV insertion events are observed four lattice positions into an island. The fifth atomic position into the step is not susceptible to insertion, so the insertion probability reaches a maxima at 15 eV. Above this energy, insertions continue to dominate beyond the first position above the step, but positions near the step become unstable and often form adatom-vacancy pairs.

Adatom-vacancy formation on the flat terrace begins abruptly near 20 eV, and the total yield increases at a rapid rate, reaching a yield of one by 60 eV. On average, more than two adatom-vacancy pairs are created per incident atom by 150 eV. By 80 eV, adatom-vacancy has a higher probability at all the step positions considered than any other mechanism.

At about 40 eV we begin to observe atoms escaping from the system, with some preference for positions close to the atomic step edge. At low energies, an atom incident just above a step could shift the registry of atoms and become incorporated, but at higher energies the transverse momentum provided by the incident atom can eject step atoms from the system. The rapid increase in yields for adatom-vacancy pairs and sputtered atoms combine to double the number of dislodged atoms between 60 an 100 eV, greatly increasing the surface mobility and reducing the net growth rate to 65% of the incident flux. More comprehensive studies of re-sputtering have been reported elsewhereWad2 .

## IV Results and Discussion of 50K KMC-MD deposition

While the molecular dynamics simulations of isolated impacts estimates the relative yields of different atomic mechanisms at various energies, the effect on a dynamically growing film requires the more sophisticated KMC-MD simulation. We have used KMC-MD to grow copper thin films on a Cu(111) surface using energies ranging from thermal to 40 eV. All the mechanisms identified with the isolated molecular dynamics in the previous section are active, but the yield for re-sputtering below 40 eV is negligible.

Five examples of copper thin films grown with the KMC-MD are presented in Figure 2. The films shown were grown with a wide range of energies: thermal (A), 12 eV (B), 21 eV (C), 30 eV (D), and 40 eV (E). All films described in this section are grown at $``$50 K at 1 monolayer/s deposition rate on an 80x80 lattice, except thermal deposition, which had a 150x150 lattice. System sizes are selected to avoid finite size effects, but all the images have the same lateral length scale, and the same color map for ease of comparison (the size of an atom is the same in all the images, and the layer depths have the same color sequence in all images). While four mono-layers of copper has been deposited in all cases, the films grown at 21 eV and 30 eV do not have any atoms in the seventh and eighth layers. The thermally deposited film has a large population of atoms in these upper two layers and is rougher than those grown with energetic deposition. The step density in the thermal films is much higher at 0.74 than the 21 eV or 30 eV films (0.39 and 0.45, respectively), corresponding to a shorter lateral length scale (step densities are discussed in more detail later).

A common way of representing smooth growth that establishes a connection with experimental effortsPoels ; Mich1  is to plot the simulated anti-phase Bragg intensities associated with refelection high energy electron diffraction or x-ray diffraction. Anti-phase intensities will exhibit complete oscillations between 0 and 1 for perfect layer-by-layer growth, and a monotonic decay for three-dimensional roughening. Due to space constraints, we have limited the presentation of simulated anti-phase intensity to our discussion of temperature and flux in Section V. With the exception of thermally deposited films, all the films studied exhibited layer-by-layer oscillations of varying strengths. The anti-phase intensity of the thermally grown film decays monotonically in this low-temperature regime, consistent with experimental observation Rosen2 . The anti-phase intensity oscillations are strongest between 20-30 eV, corresponding with the minima in roughness shown in Fig 3.

The roughness of the KMC-MD films is quantified by calculating the RMS roughness at the completion of each of the four mono-layers deposited. This roughness data as a function of the deposition energy is shown in Figure 3. Even after depositing only one monolayer, the films grown with atoms in the 20 eV range have a much lower roughness than those grown with higher or lower energies. As the film progresses, this minimum roughness appears to shift toward higher energies. The roughness of films grown with less than 25 eV grows more quickly after 1 monolayer of coverage than the roughness of the films grown with energies greater than 25 eV. Careful examination of the 40 eV data reveals very little change in the surface roughness after the first monolayer.

The time evolution of the RMS roughness is shown in greater detail for a few selected energies in Figure 4. The RMS roughness for the thermally deposited film diverges as a power law, as expectedFam . For all the energies shown, the RMS roughness grows rapidly until about 0.5 ML. Below 0.5 ML, films deposited with energies greater than 20 eV actually develop roughness faster than the thermally deposited film. As the hyperthermal beam creates large numbers of adatom vacancy pairs, the surface width increases rapidly, but these extra adatoms in turn increase nuclei densities, which contribute to higher step densities, shown in Figure 5.

Island density and step density (top frames), as well as time averaged insertion and vacancy yields (bottom frames) are presented in Figure 5A-C for 12 eV, 27 eV, and 40 eV, respectively. As the energy increases from 12 eV to 40 eV, the saturation island density (approximately the density at 0.15 ML) increases by a factor of six. Since the island density and the average island size are related ($`N_x=\theta /\overline{s}`$, where $`\theta `$ is the coverage and $`\overline{s}`$ is the average island size), one might naively expect a factor of six increase in island density to correspond with a $`\sqrt{6}`$ increase in the step density at constant coverage. In actuality, the increase in island density does not translate directly to an increase in step edge density, the step edge density in Figure 5 increases by only about 50%, not $`\sqrt{6}`$. This is partly a consequence of our definition of island density. We define the step density as the number of atoms with empty neighbor sites, since the insertion mechanism relies on displacing an atom into an empty neighboring lattice position. This definition can decouple island density and step density. For example, if all islands were composed of four atoms, the step density would be 4 times the island density. But if all the islands were made from dimers, the step density would be the same, but the island density would increase by a factor of two.

Figure 5 also displays the dynamic yields for the insertion mechanism and the formation of vacancies at 12eV (A), 27 eV (B), and 40 eV (C) respectively. While the yields discussed in the previous section and presented in Figure 1 provide an average yield for insertions and vacancies at a given energy near a (111) step edge, the yields in Figure 5 are dynamic. Each line-shape is a running average a few hundreths of a monolayer in coverage. The yields for insertions and vacancy production is observed to be sensitive to fluctuations in the surface structure.

The insertion yields at 12 eV and 27 eV track the step density very closely. At very short times, while the surface is still very flat, the step density is very low and few sites are available for insertion. During these times, the vacancy yield is large, which in turn increases the island and step density. With increasing step density, more sites become available for insertion, increasing the insertion yield. At these energies, adatom-vacancy production is suppressed at step edges. So, as the surface becomes more populated with islands and fewer flat terraces, the vacancy yield decreases. At about 0.3 ML, islands begin to coalesce and decrease the number of first layer step edges. The formation of adatom-vacancy pairs on the second layer (atop islands) keeps the step edge density high, and the insertion mechanism does not suffer.

The increase in the island density and corresponding increase in the step density below 0.5 ML in the 27 and 40 eV simulations sets the stage for smooth growth at later times. The drop in the island density to a very low value by 0.7 ML indicates near completion of the first layer before second layer growth. The abrupt change in the RMS roughness (Figure 4) near 0.5 ML of coverage illustrates the predicted benefit of using hyperthermal energy particles. With increasing particle energy the RMS roughness grows more slowly until, by 40 eV, the RMS roughness does not perceptibly increase above 1 ML of coverage. This “saturated” roughness was observed in all films grown with energies at or above 30 eV. While the roughness does not noticeably increase above 1 ML in these films, higher incidence energies result in larger “saturation” roughnesses.

The effect of a saturated roughness is the result of the insertion yield and the vacancy yield both achieving saturation. At 27 eV, the first four lattice positions immediately above a step edge are available for insertion, which suppresses vacancy formation. As a result, insertions have a higher yield during growth above 0.5 ML than vacancies (Figure 5B), and the saturation roughness is less than at 40 eV. At higher energies, the first lattice position above the step is unstable upon impact. As a result, the balance of insertion and vacancy yields falls in favor of vacancies at 40 eV (see Figure 5C, and notice the scale differences between 12 (A), 27 (B), and 40 eV (C)).

While total insertions decrease at 40 eV, other energetic effects begin to compensate and keep the roughness from increasing dramatically. For example, the increasing atom energy can break islands into smaller pieces, preventing an additional layer from nucleating on top of it. Atom impacts on top of multiple layers can lead to collective downward mobility, e.g. at 40 eV two atoms or more were observed to fall in the layer beneath the impact site one time in twenty-five. These and other mechanisms which involve collective motion of multiple atoms have been discussed in detail elsewhereVoter1 ; Wad1 .

The large insertion yield effectively reduces the interlayer diffusion barrier by providing an alternative to thermal descent for crossing the step. We have performed KMC simulations using an reduced interlayer diffusion barrier to mimic this effect. This oversimplification fails to reproduce the correct line-shape for the RMS roughness and does not provide a layer-by-layer type growth, underscoring the importance of adatom-vacancy pairs. The adatom-vacancy pairs contribute to establishing a microscopically rough interface that sustains a macroscopically smooth growth front through high insertion yields.

To review, as a film begins to grow with atoms in the 25 eV energy range, the surface initially becomes pocked with vacancies. The deposited atoms combine with adatoms from adatom-vacancy pairs to develop high densities of small islands, which have low probabilities for second layer nucleation. The vacancies and the new islands both increase the step edge density, which leads to high insertion yields. As the islands grow from insertion and aggregation, some begin to form second layer islands and obtain vacancies prior to coalescence. This establishes an average distance between step edges of about three atomic positions, and the roughness saturates. As the vacancies are filled and new levels are nucleated, the surface grows smoothly with a constant roughness. This smooth growth relies on both the insertions to keep islands growing and the vacancies to provide additional adatoms and to reduce the area available for new layer nucleation.

## V Non-intuitive Role of Temperature and Flux in Energetic Deposition

During thermal homoepitaxy, the roles of temperature and flux are well understood Rosen2 , smooth growth occurs when adatoms have enough time to diffuse to an existing island and the islands have time to coalesce before second layer nucleation occurs. This is most likely to occur when the temperature is increased to increase the diffusion length, and the flux decreased to reduce the probability of nucleating a new island before coalescence.

During hyperthermal energy deposition, this phenomenology reverses due to the strong dependence on step density. Other authors have found the best results for smooth growth with hyperthermal deposition can be obtained by maximizing the nuclei densityRosen1 ; Wulf1 . Since island densities scale as $`N_x\left(F/D\right)^p`$ (F is flux, D is the temperature dependent diffusivity, and p depends on the critical nuclei size)venables , establishing a high density of islands requires decreasing the temperature and increasing the flux. A high island density also means a small average island size ($`N_x=\theta /\overline{s}`$, where $`\theta `$ is the film coverage and $`\overline{s}`$ is the average island size), reducing the target area for second layer nucleation, and keeping the entire second layer surface close to a step edge. A hyperthermal atom incident on top one of these islands has a very high probability of inserting, rather than relying on kinetic diffusion to cross to the lower terrace.

We have presented RMS roughness as a function of time and anti-phase intensity data for films grown with 24 eV atoms at various temperatures and flux in Figure 6. Contrary to thermal depositionRosen2 , the rougher films occur with the higher temperatures and lower fluxes. In addition, one can compensate for a decrease in flux by decreasing the temperature. For example, the surface grown at 50K and 0.1 ML/s was much rougher than the film grown at 50 K and 1 ML/s, but at 35 K with 0.1 ML/s flux, the film grows smoothly.

## VI Hyperthermal Energy Induced Island Densities

While atomic insertion provides a compelling mechanism for controlling surface roughness, the extremely small island sizes and high step densities required for layer-by-layer behavior reduce the effectiveness at temperatures typical for film growth. At typical deposition temperatures, the average island size becomes large enough to reduce the number of sites available for insertion significantly, in turn reducing the average insertion yield. Even though the insertion yield drops, the adatom-vacancy pairs generated by hyperthermal energy ion beams provide additional adatoms that increase island densities.

Island densities after the first layer reaches 0.15 ML of material are shown for films grown at 50 K and 273 K in Figure 7. The 273 K deposition is performed on a 400x400 lattice at 100 ML/s using all the diffusion moves listed in Tables 1 and 2. Increased adatom density due to adatom-vacancy production leads to dramatic increases in the island density at both temperatures. The onset of sputter erosion near 40 eV acts to slow the net growth speed relative to the nucleation rate, increasing the effective flux. As material in the substrate is lost to sputter erosion, saturation island densities (typically at about 0.15 ML of coverage) are reached as low as 0.05 ML. As atoms are deposited or displaced, the first layer collects most of the material, nucleating new islands while the additional vacancies produced by the sputter erosion reduces the net meterial deposited. For example, consider an energy at which two adatom-vacancy pairs are produced and the sputter yield is 0.5. On average, every incident atom will create 2.5 atoms in the first layer, while only depositing 0.5 atoms total. After 0.1 ML of net deposition, the first layer will have almost 0.25 ML coverage.

The data presented in Figure 7 is the density of islands after the first layer has reached 0.15 ML of coverage, not 0.15 ML of total deposition. The saturation island density depends on the coverage in the layer being considered, not on the total amount of material deposited. It may be surprising such a strong effect remains at 273 K, as the thermal activation of diffusion increases the probability for adatom-vacancy recombination. We find that the high adatom densities lead to rapid formation of dimers, which still move freely with only a 117 meV diffusion barrier. Dimer step crossing is negligible, and dimer breakup occurs infrequently.

We have also found that high island densities are maintained with unexpectedly low average island sizes due to islands being “chipped” by incident atoms. While large islands are occasionally broken into two smaller stable islands, breaking dimers and adatoms off of stable islands occurs with a relatively high yield, on the order of 1/10 impacts on an island. This contributes to additional nucleation of small islands and suppresses the growth of large islands.

The significance of enhanced nuclei densities has already been experimentally demonstrated by using an ion beam to increase the free adatom density at the beginning of each monolayer of growth while depositing copper on Cu(111)Wulf1 . (In this example, a separate 1.2 keV argon ion beam and a flux of thermal atoms were used.) The ability to dramatically increase nuclei densities using a single growth beam in the room temperature regime opens new possibilities for circumventing three-dimensional growth. It is possible that pulsing the final beam energy between a low energy ($``$20 eV) and a high energy ($``$100 eV) in manner similar to the previous exampleWulf1  can allow the benefits of high island density and atomic insertion to be utilized while using a single deposition source.

## VII Acknowledgements

The authors would like to acknowledge Oana Malis, Markus Rauscher, and Joel Brock for their discussions and editorial contributions. This work was supported by the Cornell Center for Materials Research (CCMR), a Materials Research Science and Engineering Center of the National Science Foundation (DMR-0079992). Portions of this research were conducted using the resources of the Cornell Theory Center, which receives funding from Cornell University, New York State, Federal Agencies, and corporate partners. Additional support was provided by the AFOSR under grant number F49620-97-1-0020.


format: text

text:
# Lattice Predictions for Hybrids and Glueballs.

## 1 Introduction

QCD is a simple elegant theory, but it is hard to solve it using anything other than perturbation theory. Quantities, such as masses, depend on the coupling ($`g`$) like $`Me^{1/g^2}`$ , hence perturbation theory can’t be used to compute the masses of hadrons such as the proton. Our inability to solve QCD non-perturbatively makes it hard to determine basic parameters of the quark sector of the standard model. For example, the allowed range on the strange quark mass in the particle data table  is 75 to 170 MeV; a range of almost 100%.

A particularly good test of our understanding of the non-perturbative aspects of QCD is to study particles where the gauge field is excited somehow, and hence playing a more important dynamic role than in “standard” hadrons. Examples of such particles are glueballs (particles made out of the gauge fields) and hybrid mesons ($`\overline{q}q`$ and excited glue). A very good overview of the interesting issues in hadronic physics, that contrasts the hadron spectroscopy approach to the study of confinement with the results from DIS type studies, is the white paper by Capstick et al. . The only technique that offers any prospect of computing masses and matrix elements non-perturbatively, from first principles, is lattice QCD.

I review the results from the lattice for the lightest $`0^{++}`$ glueball and $`1^+`$ exotic, as these states are the closest to being experimentally confirmed. Other recent reviews  of lattice results for hybrids and glueballs focus on different aspects of the subject.

## 2 Lattice QCD calculations

Many bound state properties of QCD can be determined from the path integral

$$c(t)𝑑U𝑑\psi 𝑑\overline{\psi }\underset{\underset{¯}{x}}{}O(\underset{¯}{0},0)O(\underset{¯}{x},t)^{}e^{S_FS_G}$$
(1)
where $`S_F`$ is the fermion action (some lattice version of the continuum Dirac action) and $`S_G`$ is the pure gauge action. The path integral in eq. 1 is put on the computer using a clever finite difference formalism , due to Wilson, that maintains gauge invariance. The path integral in eq. 1 is evaluated using algorithms that are generalisations of the Monte Carlo methods used to compute low dimensional integrals. The physical picture for eq. 1 is that a hadron is created at time 0, from where it propagates to the time t, where it is destroyed. The physics from the calculation is extracted using a fit model :

$$c(t)=a_0exp(m_0t)+a_1exp(m_1t)+\mathrm{}$$
(2)
where $`m_0`$ ($`m_1`$) is the ground (first excited) state mass and the dots represent higher excitations. Although in principle excited state masses can be extracted from a multiple exponential fit, in practice this is a numerically non-trivial task, because of the noise in the data from the calculation. Any gauge invariant combination of quark fields and gauge links can be used as interpolating operators ($`O(\underset{¯}{x},t)`$) in eq. 1.

The fermion integration can be done exactly in eq. 1 to produce the fermion determinant. The determinant describes the dynamics of the sea quarks. In quenched QCD calculations, the determinant is set to a constant. Quenched calculations are roughly 1000 times cheaper computationally than the calculations that include the dynamics of the sea quarks. Quenched QCD gives quite a reasonable description of experiment. For example, the most accurate quenched calculation of the hadron spectrum, to date, has been completed by the CP-PACS collaboration . From the masses of 11 light hadrons, they conclude that the quenched approximation disagrees with experiment by at most 11%.

In an individual lattice calculation there are errors from the finite size of the lattice spacing and the finite lattice volume. State of the art lattice calculations in the quenched theory, run at a number of different lattice spacings and physical volumes and extrapolate the results to the continuum and infinite volume  limit. The increased computational costs of unquenched calculations means that most calculations are currently done at fixed lattice spacings. One of the most interesting unquenched calculations is being performed by the MILC collaboration . MILC’s calculations include 2+1 flavours of sea quarks with a lattice spacing of 0.13 fm, box size of 2.6 fm, and the lightest ratio of the pseudoscalar to vector mass is 0.4.

## 3 Results for glueballs in quenched QCD

Interpolating operators for glueballs are constructed for eq. 1 from closed loops of gauge links with the required $`J^{PC}`$ quantum numbers Some highlights of the results are that the lightest glueball is the $`0^{++}`$ state with a mass of 1.64(4) GeV . The next lightest glueball is $`2^{++}`$. The ratio of the tensor to scalar glueball mass is $`M_{2^{++}}/M_{0^{++}}`$ = 1.42(6) . The spectrum of glueball states for other $`J^{PC}`$ quantum numbers with masses under 4 GeV has been comprehensively mapped out by Morningstar and Peardon .

In the real world glueballs will decay to two mesons, hence they will have a decay width. Lattice QCD calculations are performed in Euclidean space, for convergence of the path integral in eq. 1. The Euclidean nature of lattice calculations makes the computation of inherently complex quantities such as decay widths more involved .

The GF11 lattice group computed the decay widths for the decay of the $`0^{++}`$ glueball to two pseudoscalars  to be 108(28) MeV. Although the error bar is only statistical, it is encouraging that the width was small relative to the mass, so the $`0^{++}`$ glueball may exist as a well defined state. The decay widths for individual meson pairs  did not agree with the predictions from the “flavour democratic” assumption.

The experimental situation  for light $`0^{++}`$ scalars is very interesting, because there are too many states to put into SU(3) nonets, as other particles with different $`J^{PC}`$ quantum numbers, such as the pseudoscalars, can be. The hadrons: $`f_0(1370)`$ ,$`f_0(1500)`$ , and $`f_0(1710)`$, have masses close to the mass of the quenched scalar $`0^{++}`$ glueball.

In full QCD interpolating operators with $`0^{++}`$ can be constructed from quarks and ante-quarks, such as $`\overline{\psi }\psi `$. In full QCD, the pure glue $`0^{++}`$ operators will mix with the fermionic $`0^{++}`$ operators. If the mixing is very strong, then the final $`0^{++}`$ masses will have little to do with the quenched glueball masses.

Weingarten and Lee  studied the effect of mixing between the glueball and $`\overline{\psi }\psi `$ states in quenched QCD. They measured the correlation between the $`0^{++}`$ glueball states and $`\overline{\psi }\psi `$ states in eq. 1. The results were expressed as a mixing matrix

$$\left(\begin{array}{cc}m_g& E(s)\\ E(s)& m_\sigma (s)\end{array}\right)$$
(3)
where $`m_g`$ is the glueball mass, $`m_\sigma (s)`$ is the mass of the non-singlet $`0^{++}`$ state at strange, and $`E(s)`$ is the mixing energy. Weingarten and Lee measured: $`m_g`$ = 1648(58) MeV, $`m_\sigma (s)`$ = 1322(42) MeV, and $`E(s)`$ = 61(58) MeV. The qualitative picture that emerges is that the $`f_0(1710)`$ is “mostly” $`0^{++}`$ glueball, and the $`f_0(1500)`$ is “mostly” $`\overline{s}s`$. It is not clear whether $`f_0(1500)`$ being $`\overline{s}s`$ is consistent with its decay width .

The mixing energy $`E(s)`$ has large lattice spacing errors. For example at a lattice spacing of $`a^1`$ 1.2 GeV, the Weingarten and Lee  result is $`E(s)`$ 0.36 GeV. This has been checked by another group’s result  of $`E(s)`$ 0.44 GeV.

The analysis of Weingarten and Lee  depends on the $`0^{++}`$ states being well defined in quenched QCD. Eichten  et al. have shown that there is a problem with the non-singlet $`0^{++}`$ correlator in quenched QCD. The problem can be understood using quenched chiral perturbation theory. The non-singlet $`0^{++}`$ propagator contains an intermediate state of $`\eta ^{}\pi `$. The removal of fermion loops in quenched QCD has a big effect on the $`\eta ^{}`$ propagator. The result is that a ghost state contributes to the scalar correlator, that makes the expression in eq. 2 inappropriate to extract masses from the calculation. Eichten et al.  predict that the ghost state will make the $`a_0`$ mass increase as the quark mass is reduced below a certain point. This behaviour was observed by Weingarten and Lee  for small box sizes (L $``$ 1.6 fm) for quark masses below strange. It is not clear how the problem with the non-singlet $`0^{++}`$ correlator in the quenched approximation effects the results of Weingarten and Lee , however their most important results come from masses above the strange quark mass where the ghost diagram will make a smaller contribution that may be negligible

## 4 Results for glueball masses in two flavour QCD

The Weingarten and Lee  analysis predicted that the mixing of the $`0^{++}`$ glueball and $`\overline{\psi }\psi `$ states is small. Parts of their calculation have been criticised in , however, the problems with the non-singlet $`0^{++}`$ correlator  in the quenched QCD will make further progress in mixing in the quenched QCD difficult.

A lattice QCD calculation that included the dynamics of the sea quarks should just reproduce the physical spectrum of $`0^{++}`$ states. Some insight into the composition of individual $`0^{++}`$ states, such as whether a physical particle couples to $`\overline{\psi }\psi `$ or pure glue operators, could be studied by looking at the effect of decreasing the sea quark mass. For very heavy sea quark masses the theory is more like quenched QCD, where glueballs are distinct from $`\overline{\psi }\psi `$ operators.

Hart and Teper  found that the ratio of the $`0^{++}`$ glueball mass in $`n_f=2`$ QCD to the quenched QCD result was: $`M_{n_f=2}^{0^{++}}/M_{quenched}^{0^{++}}`$ = $`0.84\pm 0.03`$ at a fixed lattice spacing of 0.1 fm. The $`n_f`$ = 2 results  for the mass of the $`2^{++}`$ were consistent with the quenched value. As the lattice spacing dependence of the mass of the singlet $`0^{++}`$ state in two flavour QCD and quenched QCD could be different, a definitive result will only come after a continuum extrapolation of the unquenched masses. In quenched QCD , the difference between the continuum extrapolated mass of $`0^{++}`$ glueball mass and the mass at 0.1 fm is of the order of 200 MeV. This is the same magnitude of the mass splittings between the masses of the experimentally observed particles $`f_0(1500)`$ and $`f_0(1710)`$. Although the current results for singlet $`0^{++}`$ states are starting to be interesting, the lattice spacing used in unquenched calculations must be reduced before direct contact can be made to phenomenology.

The mass of the $`0^{++}`$ glueball on the UKQCD data set are just degenerate with the mass of two pions . As the mass of the sea quarks is reduced, two pion states may effect the physics of singlet $`0^{++}`$ states.

## 5 Results for light $`1^+`$ exotic mesons

The quark model predicts the charge conjugation (C = $`(1)^{L+S}`$ ) and parity (P= $`(1)^{L+1}`$ ) of a meson with spin $`S`$ and orbital angular momentum $`L`$. States with quantum numbers not predicted by the quark model, such as: $`J_{exotic}^{PC}`$ = $`1^+`$, $`0^+`$, $`2^+`$, $`0^{}`$ are known as exotics . Exotic states are allowed by QCD. Morningstar and Peardon  claim that there are no glueballs with exotic quantum numbers with masses less than 4 GeV.

There are a number of different possibilities for the structure of an exotic state. An exotic signal could be: a hybrid meson, that is a quark and anti-quark with excited glue, or bound state of two quarks and two anti-quarks ($`\overline{\psi }\overline{\psi }\psi \psi `$).

One possible interpolating operator , that can be used in eq. 1, for a hybrid $`1^+`$ particle is

$$O_{1^+}(\underset{¯}{x},t)=\overline{\psi }(\underset{¯}{x},t)\gamma _jF_{ij}(\underset{¯}{x},t)\psi (\underset{¯}{x},t)$$
(4)
where $`F`$ is the QCD field strength tensor. If $`F`$ is removed from eq. 4, the operator creates the $`\rho `$ particle. In this formalism a gauge invariant interpolating operator, for any possible exotic hybrid particle or four particle state can be constructed. The dynamics then determines whether the resulting state has a narrow decay width, hence it can be detected experimentally. In the large $`N_c`$ (number of colours) limit  both exotic hybrid mesons and non-exotic mesons have widths that are small compared to their masses.

There have not been many new calculations of the mass of the light $`1^+`$ hybrid recently. All the results from the various lattice QCD calculations, by UKQCD , MILC  and SESAM  are essentially consistent with the mass of the $`1^+`$ state around $`1.9(2)`$ GeV . The interpolating operators used to create the exotic meson states in the MILC calculations  are different to those used in the UKQCD  and SESAM simulations , hence giving confidence that the systematic errors are under control. The results for the hybrid masses reported by Lacock and Schilling , include some effects from dynamical sea quarks. The recent results for the $`1^+`$ mass from calculations that used an asymmetric  lattice in time are consistent with the older results.

There are a number of experimental candidates for light $`1^+`$ states . The E852 collaboration have reported  a signal for $`1^+`$ state around $`1.6`$ GeV. There is also an experimental signal for a $`1^+`$ state at $`1.4`$ GeV .

There has been some recent work  on the quark mass dependence of the $`1^+`$ states. The lattice calculations are usually done at large quark masses and the results extrapolated to the physical quark masses. The conclusion of  was that the inclusion of the decay of the hybrid in the quark mass dependence of the exotic mass could reduce the final answer by 100 MeV. The predictions in  will be tested as the quark masses used in lattice calculations are reduced. The $`1^+`$ state at 1.4 GeV seems low relative to the lattice results.

It is possible that the states seen experimentally are really $`\overline{\psi }\overline{\psi }\psi \psi `$ states, in which case the operators used in the lattice simulations (eq. 4) might not couple strongly to them. Alford and Jaffe  studied $`\overline{\psi }\overline{\psi }\psi \psi `$ operators with $`J^{PC}`$ = $`0^{++}`$ in a recent lattice calculation. The motivation was to gain insight into states such as the $`f_0(980)`$ that some people believe is not a $`\overline{\psi }\psi `$ meson, but a $`\overline{\psi }\overline{\psi }\psi \psi `$ state. A similar lattice calculation could in principle be done for the $`J^{PC}`$ = $`1^+`$ exotic.

To definitely identify a particle requires both the calculation of the mass as well as the decay widths. There has been very little work on strong decays on the lattice. The most obvious hadronic process to study using lattice gauge theory is the $`\rho \pi \pi `$ decay, however there have only been a few attempts to calculate the $`g_{\rho \pi \pi }`$ coupling . Michael discusses the problems with the formalism for hadronic decays on the lattice .

In the static quark limit the exotic states on the lattice are described by adiabatic potentials. The ground state of the static potential ($`A_{1g}`$) is the familiar Coulomb plus linear potential. The excited potential ($`E_u`$) is a very flat potential, that can be used with Schrödinger’s equation to predict the spectrum of heavy-heavy hybrids . UKQCD  have investigated the de-excitation of the $`E_u`$ potential to the $`A_{1g}`$ potential by the emission of a light quark loop. In the real world, the decays would correspond to $`1^+\chi _b\eta `$ and $`1^+\chi _bS`$ with $`S`$ a scalar and $`\eta `$ a pseudoscalar. The decay width of $`1^+\chi _b\eta `$ and $`1^+\chi _bS`$ transitions were less than 1 MeV and around 80 MeV respectively. The various approximations in the static limit mean that these widths have no direct relevance to experiment.

The MILC collaboration  have investigated the mixing between the operator in eq. 4 and the operator ($`\pi a_1`$) eq. 5.

$$\overline{\psi }^a\gamma _5\psi ^a\overline{\psi }^b\gamma _5\gamma _i\psi ^b$$
(5)
that has the quantum numbers $`1^+`$. This type of correlator is part of the calculation required to compute the decay width of the $`1^+`$ state to $`\rho `$, and $`a_1`$. The more complicated part is to use eq. 5 in eq. 1 requires some clever numerical work.

## 6 Conclusions

The glueball spectrum from quenched QCD is essentially complete. The key issues now are to quantify the mixing of the glueball states with $`\overline{\psi }\psi `$ operators and to determine decay widths.


text:
# Low Temperature Symmetry of Pyrochlore Oxide Cd2Re2O7

(July 8, 2002)

## Abstract

We report the X-ray study for the pyrochlore oxide Cd2Re2O7. Two symmetry-lowering structural transitions were observed at $`T_{\mathrm{s1}}`$=200 K and $`T_{\mathrm{s2}}`$=120 K. The former is of the second order from the ideal cubic pyrochlore structure with space group $`Fd\overline{3}m`$ to a tetragonally distorted structure with $`I\overline{4}m2`$, while the latter is of the first order likely to another tetragonal space group $`I4_122`$. We discuss the feature of the lattice deformation.

A series of pyrochlore oxides A2B2O7 show a variety of interesting conductive and magnetic properties, depending on the substitutions of A and B sites ref1 . Geometrical frustration on the pyrochlore lattice leads to an exotic ground state such as spin ice for magnetic insulators ref2 , while for metallic systems a unique heavy carrier behavior is found in spinel LiV2O4 having a pyrochlore lattice made of V atoms ref4 . In the pyrochlore structure, the A site occupied by a rare-earth or a post-transition metal cation forms a distorted (6+2) coordination, while the B site occupied by a transition metal cation forms a distorted octahedron BO6 with an equal B-O distance. The structure comprises two interpenetrating pyrochlore lattices composed of the A or B sites (see Fig. 1).

Recently, Hanawa et al. and the others reported the first superconductor Cd2Re2O7 ($`T_\mathrm{c}`$=1 K) in the family of pyrochlore oxides ref5 ; ref6 ; ref7 . The compound exhibits two more phase transitions at $`T_{\mathrm{s1}}`$=200 K and $`T_{\mathrm{s2}}`$=120 K, where electrical and magnetic properties change dramatically ref8 ; ref9 ; ref10 . The former is of the second order, while the latter is of the first order. These transitions are accompanied by structural transitions. Three phases appearing are named as phase I ($`T`$$`>`$$`T_{\mathrm{s1}}`$), II ($`T_{\mathrm{s1}}`$$`>`$$`T`$$`>`$$`T_{\mathrm{s2}}`$), and III ($`T`$$`<`$$`T_{\mathrm{s2}}`$). Phase I crystallizes in the ideal pyrochlore structure at room temperature with space group $`Fd\overline{3}m`$. It was reported in the recent structure study that the crystal system is still cubic within the experimental resolution in the wide temperature range down to $`T`$=10 K ref9 . Moreover, taking into account the extinction rule of reflections observed below $`T_{\mathrm{s1}}`$ in the single crystal XRD experiments, a possible cubic space group of $`F\overline{4}3m`$ was suggested. However, recent Re nuclear quadrupole resonance (NQR) experiments have indicated a lack of a threefold axis below $`T_{\mathrm{s1}}`$ implying that the true symmetry is lower than cubic ref11 . In addition, very recent high-resolution XRD measurements using single crystal barely detected a small splitting of cubic Bragg peaks below $`T_{\mathrm{s1}}`$, probably due to tetragonal deformation ref12 . In this letter, we report the structural study on the low temperature phases of Cd2Re2O7 by means of X-ray diffraction and discuss possible space groups for them.

A single crystal was synthesized as described in the previous paper ref5 . We took oscillation photographs using an imaging plate type Weissenberg camera (Mac Science DIP320V) equipped with a closed-cycle helium refrigerator, which is a specialized system with very low background, in the temperature range of 10 K$``$300 K. The intensity data was collected by a four-circle diffractometer (Mac Science MXC) using a 2$`\theta `$-$`\omega `$ scan mode and a CCD area detector (Bruker SMART APEX) equipped with a helium flow type cooler in the temperature range of 85 K$``$300 K. These measurements were performed using a 21 kW rotating-anode X-ray generator with a graphite monochromated Mo-$`K`$$`\alpha `$ radiation. The cooling or heating rate was fixed to 0.5 K/min around the first-order transition at $`T_{\mathrm{s2}}`$.

First, we carefully took oscillation photographs in the whole temperature range. Neither additional reflection breaking the face-centered lattice nor peak splitting of reflections were observed. However, the forbidden reflections were found below $`T_{\mathrm{s1}}`$=200 K. Figure 2 (a) shows the temperature dependence of peak intensity for two types of reflections ($`0kl:`$ $`k+l`$$``$$`4n`$, $`00l:`$ $`l`$$``$$`4n`$), which are extinguished in $`Fd\overline{3}m`$. They increase gradually below $`T_{\mathrm{s1}}`$ with decreasing temperature, indicating that the transition is of the second order. Moreover, the slight anomaly was observed for some reflections at $`T_{\mathrm{s2}}`$=120 K.

Figure 2 (b) and (c) show the temperature dependence of peak intensity for the fundamental reflections 008 and 666 with their equivalent ones, respectively, where the intensity is normalized to the value at 220 K, just above $`T_{\mathrm{s1}}`$, to correct X-ray absorption. It is noted in each case that the equivalent reflections exhibit different intensities below $`T_{\mathrm{s1}}`$ as well as below $`T_{\mathrm{s2}}`$, indicating the reduction of cubic symmetry below $`T_{\mathrm{s1}}`$. Moreover, the sudden step for 008 and their equivalent reflections at $`T_{\mathrm{s2}}`$ suggest that the transition is of the first order. Additionally, the breaking of Friedel’s Law was found below $`T_{\mathrm{s1}}`$ as I(008)$``$I(00$`\overline{8}`$) and I(66$`\overline{6}`$)$``$I($`\overline{6}\overline{6}`$6). This implies the lack of inversion symmetry. Note that Re and Cd atoms have large anomalous scattering factors. The different intensities for the equivalent reflections and the breaking of Friedel’s Law were also detected for other many reflections.

Since the transition at $`T_{\mathrm{s1}}`$ is of the second order, we can discuss the symmetry of phase II on the basis of the group-subgroup relations of space group. Figure 3 illustrates the group-subgroup diagram from $`Fd\overline{3}m`$ to other face-centered cubic or body-centered tetragonal for the maximal translational-equivalent subgroups ref13 . Judging from the observed extinction condition, which is compatible with $`F\overline{4}3m`$ among cubic space groups, and assuming a tetragonal crystal system, one can suggest a possible space group for phase II to be $`I\overline{4}m2`$ or $`I\overline{4}`$, that is, the subgroup of $`F\overline{4}3m`$. Among them, $`I\overline{4}m2`$ is the most probable space group for phase II, because recent convergent-beam electron diffraction (CBED) study by Tsuda et al. ref14  found the existence of a mirror symmetry at $`T`$=145 K, which is missing in $`I\overline{4}`$. Note that, among symmetry operations present for $`Fd\overline{3}m`$, an inversion symmetry is lost in $`F\overline{4}3m`$, and additionally a threefold axis is lost in $`I\overline{4}m2`$. Taking account of no observation of the cubic Bragg peak splitting in our conventional-resolution measurement, the tetragonality of this phase is unusually small, and the resulting deformation of the pyrochlore lattice is minimal. In the $`I\overline{4}m2`$ model, we estimated the shift of Re atoms from the ideal coordinate (1/4,0,7/8). Assuming that the other atoms are on the ideal positions, the movement of Re atoms is estimated to be $`\mathrm{\Delta }`$$`x`$$``$+0.005 and $`\mathrm{\Delta }`$$`z`$$``$+0.005 (ca. 0.05Å), using the structure factor at 140 K. Figure 4 shows schematic tetrahedra made of Re atoms and atom positions. The volume of the tetrahedron changes alternatingly without tilting in $`I\overline{4}m2`$.

Next, we discuss the symmetry of phase III. We checked all the body-centered tetragonal lattices, and found that only the space groups shown in Fig. 3 could construct the pyrochlore structure. Since the transition at $`T_{\mathrm{s2}}`$ is of the first order, the space group at phase III seems not to be included in the subgroup of $`I\overline{4}m2`$. Therefore, possible space groups are $`I4_1/amd`$, $`I\overline{4}2d`$, $`I4_1/a`$, $`I4_122`$, $`I4_1md`$ and $`I4_1`$. Here, we can exclude $`I4_1/amd`$, $`I\overline{4}2d`$, $`I4_1/a`$ and $`I4_1md`$, because these are not compatible with the observed extinction condition of reflections with considering a possible crystal twinning. On the other hand, there is no way to distinguish a true one among the other two space groups of $`I4_122`$ and $`I4_1`$ from our experiments. However, the CBED experiments suggested that the point group of phase III is 422 ref14 . This fact makes one space group $`I4_122`$ selected unambiguously. Thus, it is plausible that a symmetry change from $`I\overline{4}m2`$ to $`I4_122`$ occurs at $`T_{\mathrm{s2}}`$, which means missing a mirror symmetry and adding a twofold axis. Such an exchange in symmetry operations is generally seen in a first-order phase transition. It is interesting to note that the Re atom coordinate changes from (x$``$1/4,0,z$``$7/8) in $`I\overline{4}m2`$ to (1/4,x$``$0,7/8) in $`I4_122`$. Then, in the $`I4_122`$ structure, the three kinds of Re-Re bond lengths exist, and the each volume of Re tetrahedon is equivalent. Let’s call the three bonds as 1, 2 and 3 in order of bond length. When we trace the bonds of each tetrahedral face as 1$``$2$``$3, as shown in Fig. 4, a left-rotation (counter-clockwise) configuration can be defined. Then, all the tetrahedral faces possess the same left-rotation configuration. It is naturally expected that a right-rotation configuration is also realized with equal probability. We call this topological chirality on the Re tetrahedra ”$`bondchirality`$” by analogy with the $`R`$$``$$`S`$ convention for a chiral molecule in organic chemistry. The direction of rotation depends on the sign of the Re atom parameter. In a real crystal, two types of domain with the left and right rotations are formed below $`T_{\mathrm{s2}}`$ as a racemic mixture. The Re tetrahedron in $`I\overline{4}m2`$ is considered to be a disordered state without bondchirality. This means that the transition at $`T_{\mathrm{s2}}`$ is relevant to ordering of the bondchirality.

In summary, we have reported the low temperature X-ray study for the first pyrochlore oxide superconductor Cd2Re2O7. The breaking of the cubic lattice and inversion symmetry are observed below $`T_{\mathrm{s1}}`$=200 K as well as below $`T_{\mathrm{s2}}`$=120 K. From the experimental findings and the group-subgroup relations, it is suggested that the space group changes from $`Fd\overline{3}m`$ to $`I\overline{4}m2`$ at $`T_{\mathrm{s1}}`$, and to $`I4_122`$ at $`T_{\mathrm{s2}}`$. In addition, the characteristic deformations of Re tetrahedron were discussed in each phase.

We are grateful to K. Tsuda, T. Kamiyama, H. Sawa, M. Takigawa, and H. Harima for helpful discussion. This work is supported by a Grant-in-Aid for Scientific Research (C) (No. 12640307) and (A) (No. 407) provided by the Ministry of Education, Science, Sports, Culture and Technology, Japan.


text:
# Hamiltonian model for multidimensional epistasis

## I Introduction

Evolution takes place via natural selection, whereby random mutations which have a salutary effect on the fitness (survival probability and/or reproductive capability) of the individual persist in the population and lead to new variants; other, neutral mutations may simply be carried along since they do not affect the well being of the individual. Deleterious mutations usually affect the organism adversely, and the accumulation of too many will reduce the fitness drastically. Each “species” actually consists of a more or less narrow distribution in the phase space of all possible genetic states, and this distribution may shift, in response to environmental pressure, along definite evolutionary routes. Smith

An interesting problem is to explain the rather fast rates at which populations seem to be able to adapt to changing environments, which suggest that the fitness does not depend in a simple additive way on the effects of each independent mutation, but that there is a nonlinear relationship, or epistasis Book , between the effects of point mutations determining the fitness function. In fact, one may surmise that evolution is not a Markovian game, but that the fitness depends upon the history of the successive mutations, in other words, it is a function of the path taken in genomic space. Kondrashov  Thus, for a mutation leading to the development of fingers to be beneficial, say, one must already have had a mutation leading to the formation of limbs.

Posed in this way, this problem seems to demand an analysis that is intrinsically dynamical. Yet, it actually lends itself to a treatment in terms of statistical equilibria, with the appropriate choice of a fitness function. In this paper we aim to provide such a fitness function, and solve the resulting model for possibly coexisting equilibrium phases indicating different species.

## II The Model

Since Eigen first introduced the quasi-species model Eigen  there has been a huge amount of both analytical and numerical work on bitstring models of genetic evolution Stauffer ; Oliveira , where the genotype of an individual is represented by a string of Boolean variables $`\sigma _i`$, $`i=1,\mathrm{}N`$. If one takes the wild type, or the initial genotype, to consist of a string of 0’s, each point mutation is indicated by flipping the bit respresenting a given gene, from 0 to 1. The number of mutations $`m`$ is then the number of 1’s on the whole string, i.e., $`m=_i\sigma _i`$. The fitness is usually taken simply to be a function, albeit nonlinear, of $`m`$.

Clearly, each $`i`$’th variable can be considered as an independent direction in phase space, so that evolution takes place in an $`N`$ dimensional space, where $`N`$ is the length of the genome. The genotype is a vertex on an $`N`$-dimensional unit hypercube, and if only single flips from 0 to 1 are allowed at a time, the path of evolution is a random walk on the edges of this hypercube. One possible way in which a vector variable can be introduced is to consider the whole vector $`𝐕\{\sigma _i\}`$ as the argument of the fitness function. Since the position of each gene on this particular string can be assigned with some arbitrariness, one may then demand that the fitness is only increased relative to the wild type if the bits that flip to 1 occur sequentially. Kondrashov . Thus, $`(0,0,\mathrm{})`$, $`(1,0,\mathrm{})`$$`(1,1,0\mathrm{})`$ are in increasing order of fitness while $`(0,1,\mathrm{})`$ is less fit.

This demands that we introduce a cost function $`H`$ which depends on the state $`𝐕`$, and we have chosen the fitness $`f`$ to decrease exponentially with this cost function, viz.,

$$fe^{\beta H}$$
(1)
where $`\beta `$ is a measure of how effective the cost function is, in affecting the fitness. The fitness function $`f`$ can be identified as the Boltzmann factor in an equilibrium statistical model with the Hamiltonian $`H`$, at constant inverse “temperature” $`\beta ^1`$. Temperature may be seen as the amount of randomness, or disorder in the system, competing with the cost function in determining the fitness. The higher the temperature, or randomness, the weaker will be the effect of the cost function in determining the state of the system. Therefore we define

$$T\beta ^1$$
(2)
as the tolerance in the system.

For the cost function we will borrow a Hamiltonian introduced by Bakk at al. Bakk , in the context of protein folding, where it is of importance that the folding events take place in a prescribed order. Thus,

$$H=\lambda J\underset{m=1}{\overset{N}{}}\underset{i=1}{\overset{m}{}}\sigma _i(1\lambda )J\underset{i=1}{\overset{N}{}}\sigma _i.$$
(3)
It can be seen that for $`\lambda =0`$, the only state which is favorable is that with all $`\sigma _i=1`$, whereas for $`0<\lambda 1`$, all states with an uninterrupted initial sequence of $`1`$s of arbitrary length $`m`$ lead to improved fitness. Here $`J`$ is a measure of the strength of the interaction between the states of each of the sites (alleles), $`\sigma _i`$. Clearly, $`\beta `$ and $`J`$ will always occur together in this model, in the product $`\beta J`$, and we may simply absorb $`J`$ into the definition of $`\beta `$.

The fitness $`f`$ may be normalized to take values between $`(0,1)`$ if we devide the expression in Eq.(1) by the sum over all states $`𝐕`$, namely,

$$Z\underset{\{\sigma _i\}}{}e^{\beta H[\{\sigma _i\}]}.$$
(4)
This sum may be performed exactly, to give,

$$Z=\frac{2^Ne^{\beta \lambda N}}{2e^{\beta \lambda }}+e^{[\lambda (N1)+1]\beta }.$$
(5)
Using this result we may compute the expectation values (average values) of the quantities $`\psi _m_{i=1}^m\sigma _i`$, which we shall call, $`\mathrm{\Psi }_m=\psi _m`$, for $`m=1,\mathrm{},N`$. Clearly, $`\mathrm{\Psi }_m`$ is the probability that in equilibrium, at least $`m`$ initial loci on the genotype have switched to 1. One finds,

$$\mathrm{\Psi }_m=\frac{1}{Z}\left[\frac{2^{Nm}e^{m\beta \lambda }e^{N\beta \lambda }}{2e^{\beta \lambda }}+e^{\beta [(N1)\lambda +1]}\right].$$
(6)
Clearly, as $`\beta 0`$, $`\mathrm{\Psi }_m(1/2)^m`$, so that it is convenient to define the order parameters

$$\mathrm{\Phi }_m=\frac{\mathrm{\Psi }_m(1/2)^m}{1(1/2)^m},$$
(7)
which take values in the interval $`(0,1)`$. In Fig. 1, we present the results of a numerical evaluation of Eq.(6), for $`\lambda =0`$, and in Fig.2, for $`\lambda =1`$, as a function of $`xT/J`$, which is the ratio of the tolerance in the system to the strength of the epistatic interactions.

We find that for $`\lambda =0`$, there is a sharp transition for large $`N`$, at $`x>x_t`$, below which the genotype is completely ordered, with all $`\sigma _i=1`$, while for $`x>x_t`$, the whole population is in the the state with all $`\sigma _i=0`$. From an inspection of Eqs.(5,7), one sees that $`x_t=(N\mathrm{ln}2)^1`$. Thus there are only two possible species in this case, with no coexistence between them. However, for $`N\mathrm{}`$, the threshold itself goes to zero. (This can be mended if the strength of the second term in Eq.(3) is chosen to be $`NJ`$ rather than $`J`$.)

For $`\lambda =1`$, it can be seen that the sharp phase transition is no longer present (the nonzero value of $`\lambda `$ has an effect similar to turning on a magnetic field in a magnetic phase transition). For large $`x`$, all the $`\mathrm{\Phi }_m`$ decay exponentially, as $`e^{m/x}`$. However, there exist effective thresholds $`x_m`$, for $`m>1`$, below which there is a nonzero probability of encountering individuals with $`m`$ initial alleles switched to 1. This signifies that at any given $`x_{m+1}<x<x_m`$, there is coexistence between $`m`$ hybrid species, with the first $`nm`$ alleles in the mutated state. The probability of encountering an individual with $`nm`$ seqentially mutated alleles is in fact precisely $`\mathrm{\Phi }_m`$. We see that $`\mathrm{\Phi }_N0`$ for $`x>x_N`$, with $`x_N1/\mathrm{ln}2`$.

To further elucidate the meaning of “tolerance,” we may compute the relative variances $`v_m`$, where

$$v_m^2(\psi _m\mathrm{\Psi }_m)^2/\mathrm{\Psi }_m^2.$$
(8)
It is trivial to note that $`\psi _m^2=\psi _m`$, so that $`v_m^2=(1\mathrm{\Psi }_m)/\mathrm{\Psi }_m`$. Then it is straightforward enough to get,

$$v_m^2=\frac{1\left(e^{\beta \lambda }/2\right)^m}{\left(e^{\beta \lambda }/2\right)^m+\left(e^{\beta \lambda }/2\right)^N\left[2e^{\beta (\lambda 1)}e^\beta 1\right]}.$$
(9)
One may see from here that the behaviour of the system is determined by the critical value of $`\beta \lambda `$ at $`\mathrm{ln}2`$, and moreover, that the variance (or the departure from the ordered phase) also depends on whether $`m`$ is small or of the order of $`N`$, as can also be seen clearly from Fig. 2.

## III Conclusion

In summary, we have presented a Hamiltonian model of multidimensional epistasis which weights only certain paths in genotype space as being favorable. The model has tunable strength ($`J`$) of interactions between different genes, which can be absorbed into an overall parameter ($`\beta `$) which determines how tolerant the environment is to deviations from the wildtype, as well as a parameter which decides whether coexistence between hybrid individuals will or will not be allowed. The model exhibits a transition between two pure types as a function of $`\beta `$ for $`\lambda =0`$. For $`\lambda 0`$, low tolerances $`T=\beta ^1`$ give rise to the appearance of hybrid types, in case a given series of mutations increases the fitness.

Acknowledgements

We are grateful to Güneş Söyler for many useful discussions. AE acknowledge partial support from the Turkish Academy of Sciences.

Figure Captions

Fig. 1. The order parameters $`\mathrm{\Phi }_m`$ for $`\lambda =0`$ all differ from zero at the same transition point. Here the length of the genome is 100. There are no hybrid species.

Fig. 2. The order parameters $`\mathrm{\Phi }_m`$ for $`\lambda =1`$, with $`N=100`$. There is a set of $`N`$ effective thresholds, below which hybrid species arise, with $`m`$ sequentially mutated alleles.


text:
# 1 Introduction

## 1 Introduction

Evidence for a tiny but non-zero cosmological vacuum energy density has steadily mounted in recent years . With information from a variety of observations put together, it appears that we are living in a flat universe, such that $`\mathrm{\Omega }_\mathrm{\Lambda }+\mathrm{\Omega }_m=1`$, where $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and $`\mathrm{\Omega }_m=0.3`$ denote the density parameters associated with the vacuum energy and matter respectively. These observations pose at least two formidable theoretical challenges, namely, what physics determines the magnitude of the vacuum energy density to be of order

$$\rho _V10^{120}M_{\mathrm{Pl}}^4(10^3\mathrm{eV})^4,$$
(1.1)
and why are $`\rho _V`$ and $`\rho _m`$ in such close proximity in magnitude, considering that their evolution, at least for constant cosmological constant, can be so different. Here $`M_{\mathrm{Pl}}=2.410^{18}`$ GeV denotes the reduced Planck mass.

A certain amount of theoretical prejudice suggests that the true vacuum energy density of the universe could be exactly zero. Needless to say, how the zero value occurs remains a complete mystery, despite many attempts . Nonetheless, following this reasoning, it has been speculated that the present universe happens to lie in a false vacuum, separated from the true vacuum by the observed energy density $`\rho _V`$ . The modest goal of this approach is to try to identify the physics underlying the origin of $`\rho _V`$, as well as ensure that the metastable vacuum is sufficiently long lived.

In this letter we propose a new mechanism for realizing $`\rho _V`$ of the desired magnitude which is based on warped geometry. While the original motivation of warped geometry was to resolve the notorious gauge hierarchy problem , it has subsequently been exploited to provide at least a qualitatively new understanding of fermion mass hierarchies and their mixings, especially within the context of neutrino oscillations . We will show here that the warped geometry can be further exploited to provide a new mechanism for generating $`\rho _V`$ of the desired magnitude. Earlier attempt in this direction have been made in ref. . Recall that the TeV scale in warped geometry arises from the presence of the warp factor $`\mathrm{\Omega }^1=\mathrm{exp}(\pi kR)=\mathrm{TeV}/M_P`$, where $`1/k`$ denotes the AdS curvature radius and $`R`$ is the radius of the orbifold $`S^1/Z_2`$. We will see that $`\rho _V`$ can be associated with the vacuum energy density of a suitable scalar field, such that

$$\rho _V^{1/4}\mathrm{\Omega }^2M_{\mathrm{Pl}}=\mathrm{TeV}^2/M_{\mathrm{Pl}}10^3\mathrm{eV}.$$

We also discuss how the quintessence scenario  can be realized within this approach.

## 2 Scalar Fields in a Warped Background

We take the fifth dimension to be an $`S_1/Z_2`$ orbifold with a negative bulk cosmological constant, bordered by two 3-branes with opposite tensions and separated by distance $`R`$. Einstein’s equations are satisfied by the non-factorizable metric

$$ds^2=e^{2\sigma (y)}\eta _{\mu \nu }dx^\mu dx^\nu +dy^2,\sigma (y)=k|y|$$
(2.2)
which describes a slice of AdS5. The 4-dimensional metric is $`\eta _{\mu \nu }=\mathrm{diag}(1,1,1,1)`$, $`k`$ is the AdS curvature related to the bulk cosmological constant and brane tensions, and $`y`$ denotes the fifth coordinate. The AdS curvature and the 5d Planck mass $`M_5`$ are both assumed to be of order $`M_{\mathrm{Pl}}`$. The AdS warp factor $`\mathrm{\Omega }=e^{\pi ky}`$ generates an exponential hierarchy of energy scales. If the brane separation is $`kR11`$, the scale at the negative tension brane, located at $`y=\pi R`$, is of TeV-size, while the scale at the brane at $`y=0`$ is of order $`M_{\mathrm{Pl}}`$. At the TeV-brane gravity is weak because the zero mode corresponding to the 4D graviton is localized at the positive tension brane (Planck-brane).

Let us consider the equation of motion of a real scalar field in the background (2.2)

$$\frac{1}{\sqrt{g}}_M(\sqrt{g}g^{MN}_N\mathrm{\Phi })M^2\mathrm{\Phi }=0$$
(2.3)
which follows from the action

$$S_5=d^4x𝑑y\sqrt{g}\left(\frac{1}{2}(_M\mathrm{\Phi })(^M\mathrm{\Phi })+\frac{1}{2}M^2\mathrm{\Phi }^2\right).$$
(2.4)
Here $`g_{MN}`$ denotes the 5D metric and $`g=e^{4\sigma }`$ its determinant. In general, the 5D scalar mass consists of bulk and brane contributions

$$M^2(y)=b^2k^2+a^2k\delta (y\pi R)+\stackrel{~}{a}^2k\delta (y),$$
(2.5)
where $`b`$, $`a`$ and $`\stackrel{~}{a}`$ are dimensionless parameters. In the case of supersymmetry the relations

$`\stackrel{~}{a}^2`$ $`=`$ $`a^2`$
$`b^2`$ $`=`$ $`{\displaystyle \frac{1}{4}}a^2(a^2+8)`$ (2.6)
have to be satisfied , and $`\mathrm{\Phi }`$ becomes a complex scalar field which, together with a second complex scalar and a Dirac fermion, forms a hypermultiplet. The connection to the 5D mass parameter of the Dirac fermion $`M_\mathrm{\Psi }=c\sigma ^{}`$ is given by $`c=3/2+a^2/2`$. A special case of eqs. (2) is the massless scalar field where $`a^2=\stackrel{~}{a}^2=b=0`$ (and $`c=3/2`$).

The effective 4D theory arises from the Kaluza-Klein (KK) decomposition

$$\mathrm{\Phi }(x^\mu ,y)=\frac{1}{\sqrt{2\pi R}}\underset{n=0}{\overset{\mathrm{}}{}}\mathrm{\Phi }^{(n)}(x^\mu )f_n(y),$$
(2.7)
where the wave functions $`f_n(y)`$ satisfy the differential equation

$$(_5^2+4\sigma ^{}_5+M^2)f_n=e^{2\sigma }m_n^2f_n,$$
(2.8)
and $`m_n`$ are the masses of the KK excitations. The solution to eq. (2.8) is given by

$$f_n(y)=\frac{e^{2\sigma }}{N_n}\left[J_\alpha (\frac{m_n}{k}e^\sigma )+\beta _\alpha (m_n)Y_\alpha (\frac{m_n}{k}e^\sigma )\right],$$
(2.9)
with $`\alpha =\sqrt{4+b^2}`$. Assuming $`\mathrm{\Phi }`$ is even under the $`Z_2`$ orbifold transformation, i.e. $`f_n(y)=f_n(y)`$, the coefficient $`\beta `$ and the KK spectrum follow from the matching conditions at the branes,

$`\beta _\alpha (x_n,\stackrel{~}{a}^2)`$ $`=`$ $`{\displaystyle \frac{(\stackrel{~}{a}^2+42\alpha )J_\alpha (x_n)+2x_nJ_{\alpha 1}(x_n)}{(\stackrel{~}{a}^2+42\alpha )Y_\alpha (x_n)+2x_nY_{\alpha 1}(x_n)}},`$ (2.10)
$`\beta _\alpha (x_n,\stackrel{~}{a}^2)`$ $`=`$ $`\beta _\alpha (\mathrm{\Omega }x_n,a^2),`$ (2.11)
where $`x_n=m_n\mathrm{\Omega }/k`$. Note that for non-vanishing boundary mass terms the derivative of $`f_n`$ becomes discontinuous on the boundaries. The normalization constants $`N_n`$ in (2.9) are defined such that

$$\frac{1}{2\pi R}_{\pi R}^{\pi R}𝑑ye^{2\sigma (y)}f_n^2(y)=1.$$
(2.12)
The mass splitting between the KK modes is $`\mathrm{\Delta }m_n\pi k\mathrm{\Omega }^1`$ TeV.

In the massless case the zero mode wave function is constant, i.e. $`f_0(y)=1/N_0\sqrt{2\pi kR}`$. The KK states have TeV scale masses and are localized at the TeV-brane. The mass of the first excited state is $`m_13.83k\mathrm{\Omega }^1`$. The effects of the bulk and brane mass terms on the mass and wave function of the lowest mode are quite different. If one switches on a small bulk mass, one finds $`m_0bk/\sqrt{2}`$ for $`b\stackrel{<}{}\mathrm{\Omega }^1`$ ($`bk\stackrel{<}{}`$ TeV). Thus, there is no warp factor suppression. The zero mode wave function is still approximately constant, with a small enhancement at the TeV-brane. Increasing the bulk mass to $`b\stackrel{<}{}1`$ ($`bk\stackrel{<}{}M_{\mathrm{Pl}}`$) the mass of the lightest mode remains frozen at $`x_0=3.83`$, which is practically the former first excited state. For $`b\stackrel{>}{}1`$ ($`bk\stackrel{>}{}M_{\mathrm{Pl}}`$) one finds that $`m_0`$ also increases. The mass term at the Planck-brane induces a similar behavior of the zero mode, except there is no increase of $`m_0`$ for $`\stackrel{~}{a}\stackrel{>}{}1`$ ($`\stackrel{~}{a}k\stackrel{>}{}M_{\mathrm{Pl}}`$).

In the case of a TeV-brane mass the result is very different. Because of the factor $`\sqrt{g}`$ in the action (2.4), the mass of the zero mode is warped by two powers of $`\mathrm{\Omega }`$. With this double warping we have

$$m_0ak\mathrm{\Omega }^2,a\stackrel{<}{}1.$$
(2.13)
The wave function is nearly constant, with a dip at the TeV-brane, as shown in fig. 1a for $`kR=11.13`$ ($`M_{\mathrm{Pl}}e^{\pi kR}=1.6`$ TeV). For $`a\stackrel{>}{}1`$ the dip in the wave function slows down the increase of $`m_0`$. For large values of $`a`$ any further increase is completely compensated by this, so that the zero mode mass saturates at $`m_0=2\sqrt{2}k\mathrm{\Omega }^2`$. The dependence of $`m_0`$ on the TeV-brane mass is displayed in fig. 1b. We note that gauge bosons show a similar behavior as we discussed in ref. . Numerically, for a TeV-brane mass of order $`M_{\mathrm{Pl}}`$, $`m_0`$ is comparable to the energy scale associated with dark energy (1.1). Of course, at this stage $`m_0`$ is not yet related to some vacuum energy, but merely represents a tiny particle mass. In the next section we will show how such a link can be established.

At the end of this section, let us briefly summarize the supersymmetric case which was extensively discussed in refs. . Supersymmetry requires the bulk and brane mass terms to be simultaneously present. Because of the special relations eq. (2), the zero mode remains strictly massless, despite the non-zero bulk and brane masses. The wave function has an exponential form ,

$$f_0(y)=\frac{1}{N_0}e^{\frac{a^2}{2}\sigma }.$$
(2.14)
For a positive (negative) mass squared term at the TeV-brane the scalar field is localized towards the Planck-brane (TeV-brane). The supersymmetric scalar fields thus behave similar to fermions in the warped background, which also possess massless exponential zero modes in the presence of bulk mass terms .

## 3 Exponentially Small Contributions to the Scalar Potential

In the last section we learned that a Planck-size mass term at the TeV-brane is translated into a sub-eV mass of the scalar particle. We now generalize this result to include a scalar potential on the TeV-brane. This will allow as to generate tiny energy splittings in the minima of the 4D effective potential, which we hope to identify with the observed dark energy.

We start with the simplest case of a massless scalar field with $`a=b=\stackrel{~}{a}=0`$ which, for instance, could be a consequence of supersymmetry. If we now switch on some general potential on the TeV-brane, related for example, to supersymmetry breaking, the scalar zero mode remains, as discussed in the last section. Consider the following contribution to the 5D action

$$\mathrm{\Delta }S_5=d^4x𝑑y\sqrt{g}\left(\frac{1}{2}a^2k\mathrm{\Phi }^2+\frac{1}{3}E_5k^{1/2}\mathrm{\Phi }^3+\frac{1}{4}\lambda _5k^2\mathrm{\Phi }^4\right)\delta (y\pi R).$$
(3.15)
We approximate $`f_0`$ by a constant, and after integrating over the extra dimension, we arrive at the 4D effective potential for the zero mode $`\varphi `$

$$V_4=\frac{1}{2}m^2\varphi ^2+\frac{1}{3}E\varphi ^3+\frac{1}{4}\lambda \varphi ^4,$$
(3.16)
where

$`m^2`$ $`=`$ $`a^2k^2\mathrm{\Omega }^4`$
$`E`$ $`=`$ $`E_5k\mathrm{\Omega }^4`$
$`\lambda `$ $`=`$ $`\lambda _5\mathrm{\Omega }^4.`$ (3.17)
Thus all operators on the TeV-brane are red-shifted by four powers of the warp factor. The warped geometry not only generates tiny masses but also exponentially small coupling constants. This result may seem surprising since the zero mode $`f_0`$ is constant in the extra dimension and thus has a large overlap with the TeV-brane. However, as we observe from eq. (2.12), in the integration over the extra dimension the region close to the TeV-brane is exponentially suppressed. The rescaled zero mode $`\stackrel{~}{f}_0=e^\sigma f_0`$ is the more appropriate object to look at. For a constant $`f_0`$, the rescaled zero mode is exponentially localized towards the Planck-brane, which explains the tiny quantities in eq. (3). For a very large mass squared term on the TeV-brane, i.e. $`a1`$, there is an additional suppression of the brane operators (3.15) induced by the dip in the zero mode. We note that scalar couplings in the bulk or on the Planck-brane are not scaled down by the warped geometry.

We now assume that $`m^2<0`$ in order to generate a vev for the scalar field. Without additional interactions from the bulk, the tiny self coupling $`\lambda `$ induced by the brane can stabilize the field only at a Planck-size vev $`\varphi k`$. We therefore have to introduce some additional interaction. A quartic coupling in the bulk does not feel the warping, and therefore generates $`\lambda 1`$ in the effective 4D action. A quartic coupling on the Planck-brane would have the same effect. The scalar field then acquires a tiny vev $`\varphi k\mathrm{\Omega }^2`$. If there is now more than one minimum with a similar vev, the energy splitting between the two minima will be on the order of $`\mathrm{\Delta }Vm^2\varphi ^2k^4\mathrm{\Omega }^8`$. This is exactly what is required for the dark energy. Two minima can easily arise, for example, in the case of two interacting scalar fields with a polynomial potential. We will show later that the false vacuum has a lifetime much longer than the age of the universe. In the next section we will also discuss how this setup can be stabilized against radiative corrections in a supersymmetric framework.

A different possibility for the scalar field to be stabilized is by some non-renormalizable interaction in the bulk, for instance, by a $`(1/Q^2)\varphi ^6`$ operator. As discussed in refs. , non-renormalizable interaction may be suppressed by much smaller energy scales than $`M_{\mathrm{Pl}}`$ because of the warped geometry. However, as explained above, in some sense the flat zero mode is localized towards the Planck-brane. Therefore non-renormalizable operators are not enhanced, and we expect $`QM_{\mathrm{Pl}}=k`$ in the 4D action. The scalar would then be stabilized at $`\varphi =\sqrt{mQ}k\mathrm{\Omega }^1`$, which is the KK scale. Assuming that the $`Z_2`$ symmetry in the scalar potential is broken by the cubic term on the brane, the induced energy splitting between the vacua at $`\pm \sqrt{mQ}`$ is $`\mathrm{\Delta }VE\varphi ^3k^4\mathrm{\Omega }^7`$. This is a factor of $`\mathrm{\Omega }`$ larger than what we found in the previous example. To make it compatible with the observational value (1.1) we either would have to increase the warp factor and end up with a KK scale of about 10 GeV(!), or to reduce the suppression scale to about $`Q=k\mathrm{\Omega }^{1/3}`$. Both possibilities do not look realistic. In the next section we discuss the supersymmetric case (with a non-vanishing bulk mass) which gives us additional possibilities.

The false vacuum energy not only leads to an exponential expansion of the universe, but will also modify the geometry of the extra dimension. To estimate the size of the back reaction we compare the false vacuum energy density to the bulk cosmological constant and brane tensions of the unperturbed Randall-Sundrum solution . We consider a somewhat simplified setup where the scalar vev is induced by a negative squared mass term $`(1/2)a^2k\mathrm{\Phi }^2`$ at the TeV-brane, as in eq. (3.15). The potential is stabilized by a quartic self interaction in the bulk, $`(1/4)\lambda _\mathrm{b}\mathrm{\Phi }^4`$. Because the scalar potential changes along the extra dimension, the vacuum configuration $`\mathrm{\Phi }(y)`$ does as well and has to be computed numerically. We find the profile to be $`\mathrm{\Phi }(y)ak\mathrm{\Omega }^2/\sqrt{\lambda _\mathrm{b}}k^{3/2}\mathrm{\Omega }^2`$ and almost flat in the extra dimension. In the second step we assumed $`a1`$ and $`\lambda _\mathrm{b}k^1`$. Close to the TeV-brane, where the instability occurs, the profile $`\mathrm{\Phi }(y)`$ increases by about 10 percent relative to its almost constant value in the bulk. On the TeV-brane the scalar potential is $`V=(1/2)a^2k\mathrm{\Phi }^2k^4\mathrm{\Omega }^4`$. In the bulk we obtain $`V=(1/4)\lambda _\mathrm{b}\mathrm{\Phi }^4k^5\mathrm{\Omega }^8`$. This means that Einstein’s equations are still completely dominated by Planck-size bulk cosmological constant and brane tensions of the Randall-Sundrum solution . Thus we can safely ignore the back reaction, at most of order $`\mathrm{\Omega }^4`$, on the wave functions (2.9) and (2.14).

## 4 The Supersymmetric Case

A crucial point is that a perturbation of the mass relations (2) by some supersymmetry breaking mass term on the TeV-brane with $`\mathrm{\Delta }a^2\stackrel{<}{}1`$ leaves the zero mode (2.14) practically unchanged. For $`\mathrm{\Delta }a^2\stackrel{>}{}1`$ we observe an additional reduction of $`f_0`$ close to the TeV-brane, as in the case of a vanishing bulk mass. To obtain the 4D couplings from the brane action (3.15) we approximate the zero mode by its supersymmetric shape (2.14). We find

$`m^2`$ $``$ $`\mathrm{\Delta }a^2k^2\mathrm{\Omega }^{4a^2}`$
$`E`$ $``$ $`E_5k\mathrm{\Omega }^{4\frac{3}{2}a^2}`$
$`\lambda `$ $``$ $`\lambda _5\mathrm{\Omega }^{42a^2},`$ (4.18)
where $`a^2`$ is now the TeV-brane mass term required by supersymmetry and determines the shape of the zero mode. We have omitted a factor $`N_0\sqrt{2\pi kR}`$ to some power which is of order unity. Localizing the scalar zero mode towards the Planck-brane by taking $`a^2>0`$, we can suppress the induced 4D couplings even stronger than in eq. (3). This opens up new possibilities to generate a small vacuum energy density.

The ($`N=2`$) 5D supersymmetry prohibits a self interaction of the scalar field in the bulk. Self interactions can be provided, however, by superpotentials on the branes, where the boundary conditions break the supersymmetry to $`N=1`$. Planck-brane operators are not warped down and induce order unity couplings in the 4D potential, in contrast to their TeV-brane counterparts (4).

Localizing the scalar field close to the Planck-brane, we can implement a quintessence-like solution to the dark energy problem. The cosmological evolution is assumed to leave the scalar field displaced from its minimum towards which it is slowly evolving. For polynomial potentials, the necessary slow rolling conditions are satisfied only for very large vevs $`\varphi M_{\mathrm{Pl}}`$. A small vacuum energy can therefore be explained only by extraordinarily small couplings. To generate $`\rho _V`$ from a term $`m^2\varphi ^2m^2M_{\mathrm{Pl}}^2`$ requires $`m^2\mathrm{\Omega }^8M_{\mathrm{Pl}}^2`$. From eq. (4) we see that this requires $`a^24`$. We can also implement this scenario by relying on the quartic term of the brane potential, $`\rho _v\lambda \varphi ^4`$. For $`a^22`$ we obtain the required small $`\lambda \mathrm{\Omega }^8`$. The other terms in eq. (3.15) as well as possible bulk self interactions have to vanish to accommodate this solution. It remains to be seen whether a realistic scenario of this type can be constructed in which the required small masses remain stable under quantum corrections .

A true cosmological constant, instead of quintessence, can be generated from a supersymmetry breaking cubic coupling which shifts previously degenerate minima. We assume that a scalar self interaction on the Planck-brane leads to $`\lambda 1`$ in the 4D action. If the scalar potential is now destabilized by a supersymmetry breaking brane mass (4), the induced vev is $`\varphi mk\mathrm{\Omega }^{2a^2/2}`$. The degeneracy of the minima at $`\pm \varphi `$ is lifted by the supersymmetry breaking cubic term on the brane (4). The induced energy splitting between the vacua is then of order $`\mathrm{\Delta }VE\varphi ^3k^4\mathrm{\Omega }^{103a^2}`$. The preferred value of the cosmological constant (1.1) then leads to $`a^22/3`$, i.e. the scalar field is localized somewhat towards the TeV-brane. The mass of the scalar particles associated with $`\mathrm{\Phi }`$ is of order $`\varphi k\mathrm{\Omega }^{5/3}10^2`$ eV.

## 5 Quantum corrections

The case with a vanishing bulk mass discussed in section 3, can readily be embedded in the supersymmetric framework by setting $`a=0`$ in eq. (4). Then supersymmetry helps to tame radiative corrections. As in four dimensions the soft mass of the scalar zero mode, $`m^2`$, will receive a 1-loop quantum correction $`\delta m^2(1/16\pi ^2)m^2\lambda \mathrm{ln}(\mathrm{\Lambda })`$ by the exchange of the scalar zero mode and its superpartner. Here $`\mathrm{\Lambda }`$ denotes the momentum cut-off and $`\lambda `$ stems from the self coupling localized at Planck-brane. As discussed above, in the case $`a=0`$ we have $`m^2k^2\mathrm{\Omega }^4`$ and $`\lambda 1`$. In the warped model there are additional radiative corrections to $`m^2`$ by the exchange of KK states of the scalar field. The KK states are localized towards the TeV-brane. Therefore they acquire a larger soft mass of order $`k^2\mathrm{\Omega }^2`$. The quartic coupling between two zero modes and two KK states of the scalar field is of order $`\mathrm{\Omega }^2`$. Thus loops with KK states of the scalar are of the same order as zero mode loop. It remains to be seen if performing the sum over the KK contributions, which are individually small, leads to a destabilization of the soft mass of the zero mode.

Since we assume the scalar field to be a gauge singlet, further quantum corrections can only come from gravity. Exchange of a zero mode gravitino in the loop gives rise to $`\delta m^2(1/16\pi ^2)m_{3/2}^2(\mathrm{\Lambda }/M_{\mathrm{Pl}})^2`$, where $`m_{3/2}`$ is the gravitino mass. Since for the zero mode of the gravitino $`m_{3/2}k\mathrm{\Omega }^3`$, we obtain a tiny correction even for a Planck-size cut-off. For the KK gravitinos the situation is different. Being localized towards the TeV-brane, their supersymmetry breaking mass is TeV-size $`(k\mathrm{\Omega }^1)`$  and their dimension-six interaction with the scalar zero mode is suppressed by $`1/(\mathrm{\Omega }M_{\mathrm{Pl}}^2)`$. The corresponding correction to the scalar mass is then $`\delta m^2(1/16\pi ^2)k^2\mathrm{\Omega }^3(\mathrm{\Lambda }/M_{\mathrm{Pl}})^2`$. Since the dimension-six operator has its support at the TeV-brane, it seems plausible that the KK gravitino loop is cut off at the TeV-scale. In this case its contribution to the scalar mass $`\delta m^2(1/16\pi ^2)k^2\mathrm{\Omega }^5`$ is safely suppressed. However, given the limitations and uncertainties which are inherent in our treatment of radiative corrections, the issue of the quantum stability of our approach to the dark energy problem is not yet satisfactorily settled.

## 6 Lifetime of the False Vacuum

The false vacuum state is metastable and will finally relax to the true vacuum by thermal fluctuations and/or quantum tunneling. We must require the lifetime of the false vacuum to be longer than the age of the universe. Because of the low temperature of the present universe, quantum tunneling is the dominant decay process. In the WKB approximation the tunneling probability is given by

$$p\left(\frac{t_U}{R}\right)^4e^B,$$
(6.19)
where $`t_U`$ denotes the age of the universe, and $`R1/m`$ is the characteristic scale of the problem. In the thin wall approximation the Euclidean action of the tunneling configuration

$$B=27\pi ^2\frac{S_1^4}{2(\mathrm{\Delta }V)^3},$$
(6.20)
depends on the energy splitting of the vacua $`\mathrm{\Delta }V`$ and the “surface tension” of the bubble $`S_1=𝑑\varphi \sqrt{2V(\varphi )}`$.

Vacuum decay is most effective between minima which are close to each other in field space and separated only by a low energy barrier $`V_B`$. Among the examples discussed before the tunneling rate is largest in the case of section 3 where $`m\varphi 10^3`$ eV. From the constraint $`p<1`$ we obtain $`B\stackrel{>}{}280`$. For a polynomial potential one finds $`B710^3(V_B/\mathrm{\Delta }V)^3/\lambda `$, i.e. the false vacuum is metastable even for $`V_B\mathrm{\Delta }V`$, as we expect for a potential with only a single mass scale. In the other example we discussed, the barrier heights and scalar vevs are even much larger. Thus, in all the cases discussed the false vacuum is sufficiently long lived.

## 7 Conclusions

In conclusion, we have shown that a warped geometry setting could help provide a new explanation of the origin of the tiny vacuum energy density of the universe that is indicated by recent observations. This can either be in the form of a cosmological constant or, if supersymmetry is invoked, as contribution arising from a slow rolling quintessence field. Although no explanation exists as to why all other contributions to the vacuum energy density effectively vanish, it is certainly intriguing to think that resolution of the gauge hierarchy and dark energy problems may have a common origin.

## Acknowledgements

We acknowledge useful discussions with A. Hebecker and A. Pomarol and a helpful correspondence with T. Gherghetta. We wish to thank the Alexander von Humboldt Stiftung for providing the impetus for our collaboration. Q.S. also acknowledges the hospitality of the Theory Groups at DESY and Heidelberg, especially Wilfried Buchmüller, Michael Schmidt and Christof Wetterich. This work was supported in part by DOE under contract DE-FG02-91ER40626.