text:
# Identifying structural patterns in ‘disordered’ metal clusters

## I Introduction

There has been much recent theoretical interest in the possibility that small clusters could have lowest-energy structures that are disordered or amorphous. Examples have been found for model clusters interacting with long-ranged pair potentials,Doye et al. (1995); Doye and Wales (1997) and a wide variety of metal clusters.Glossman et al. (1993); Garzón et al. (1998); Doye and Wales (1998); Garzón et al. (1999); Michaelian et al. (1999); Garzón et al. (2000); Wilson and Johnston (2000); Taneda et al. (2001); Michaelian et al. (2002); Oviedo and Palmer (2002); Massen et al. (2002); Darby et al. (2002); Wang et al. (2002); Doye and Hendy (2003) For some of these examples, these disordered clusters appear at sizes between the magic numbers for the usual icosahedral, decahedral or face-centred-cubic (fcc) forms, that are typically most stable for materials that are close-packed for bulk;Doye and Wales (1998); Wilson and Johnston (2000); Darby et al. (2002); Massen et al. (2002) i.e. the disordered structures are more stable than ordered structures with incomplete outer layers. More interesting are the more extreme examples, where disordered structures are lowest in energy even at the magic numbers of the common structural forms.Garzón et al. (1998, 1999); Michaelian et al. (1999, 2002) These results would suggest for these clusters that the disordered structures are dominant for the relevant size ranges.

By disordered it is usually meant that the structure has no discernible overall order. Of course, there is still local order, as evidenced by structural probes, such as the radial distribution function, but these usually have forms similar to that for bulk liquids and glasses. Additionally, another common features for these clusters is that there are many other disordered structures that lie very close in energy to the global minimum.

Disordered clusters are most likely to occur when the constraints on the nearest-neighbour distances are weak. Then, the energetic cost for the strain present in the individual nearest-neighbour bonds can be low enough to be offset by other advantageous features of the disordered clusters, such as a low surface energy. For the Morse clusters, this occurs when the potential is long-ranged and has a wide, soft well.Doye et al. (1995); Doye and Wales (1997) For metal clusters, the many-body part of the potential is relatively insensitive to disorder in the nearest-neighbour distances.Soler et al. (2000) If the contraction at the surface of the cluster (again due to the many-body forces) is strong, the ordered structures can be sufficiently destabilized to make the disordered structures lowest in energy.Soler et al. (2000)

For most of the examples where disordered structures have been found to be lower in energy than the most stable icosahedral, decahedral and fcc clusters, only a few sizes have usually been considered.Garzón et al. (1998, 1999); Michaelian et al. (1999, 2002) So, the possibility remains that at some other sizes the local structural preferences present in the disordered clusters could be assembled in a way that gives rise to a particularly stable structure with overall order. If this were the case, global optimization over a complete range of sizes would then reveal these new structural motifs and magic numbers. This has recently been done for model lead clusters, where, although none of the usual forms are ever lowest in energy and many of the clusters have no apparent overall order, particularly stable high-symmetry clusters were found at some sizes.Doye and Hendy (2003)

Indeed there are general grounds for expecting high symmetry structures to be more prevalent among the global minima, because such structures are likely to have more extreme values of the energy (both high and low).Wales (1998) This expectation, in my experience at least, seems to be born out empirically. It is rare that the global optimization of clusters does not reveal ordered high symmetry forms at some sizes. For example, even for potentials which have been designed to favour glassy configurations, global optimization has revealed the presence of unusual, but nevertheless ordered, structures.Doye et al. (2001, 2003)

Here, I wish to examine some prototypical metal clusters that so far have only been found to exhibit disordered structures. For potentials of the Gupta form the dependence of the tendency to disorder on some of the parameters of the potential has been elaborated.Michaelian et al. (1999); Soler et al. (2000) Of the parameterized metals zinc and cadmium clusters emerged as those with the strongest preference for disorder. This tendency was born out in global optimization studies at sizes where particularly stable fcc, decahedral and icosahedral structures were possible; as expected, the resulting structures appeared disordered.Michaelian et al. (2002)

In this paper, I will attempt to identify structural patterns for these two cluster systems, in particular searching for new magic numbers and novel types of order. This aim is in a similar spirit to Ref. Soler et al., 2001, where possible structural patterns for gold clusters modelled by the Gupta potential have been suggested. To achieve this I have performed global optimization for all clusters with up to 125 atoms. A further aim is to then relate back the identified structural patterns to the form of the potential.

## II Methods

### II.1 Potential

To model the zinc and cadmium clusters I use a Gupta potentialGupta (1981) fitted by Cleri and Rosato.Cleri and Rosato (1993) The potential energy is given by

$`E`$ $`=`$ $`E_{\mathrm{pair}}+E_{\mathrm{embed}}`$ (1)
$`=`$ $`{\displaystyle \underset{i<j}{}}\varphi \left(r_{ij}\right)+{\displaystyle \underset{i}{}}F\left(\overline{\rho }_i\right),`$ (2)
where $`\varphi (r)`$ is a short-ranged pair potential, $`U(\overline{\rho })`$ is a many-body embedding (or glue) function and $`\overline{\rho }_i`$ is defined as

$$\overline{\rho }_i=\underset{j}{}\rho \left(r_{ij}\right),$$
(3)
where $`\rho (r)`$ is an “atomic density” function.

For potentials of the Gupta form

$`\varphi \left(r\right)`$ $`=`$ $`2Ae^{p(r/r_01)}`$ (4)
$`F\left(\overline{\rho }\right)`$ $`=`$ $`\xi \sqrt{\overline{\rho }}`$ (5)
$`\rho \left(r\right)`$ $`=`$ $`e^{2q(r/r_01)}.`$ (6)
These forms arise from the second moment approximation of a tight-binding Hamiltonian. However, these functions are non-unique. Functions that give exactly the same energy can be constructed by the transformation

$`\varphi ^{}\left(r\right)`$ $`=`$ $`\varphi \left(r\right)+2g\rho (r)`$ (7)
$`F^{}\left(\overline{\rho }\right)`$ $`=`$ $`F\left(\overline{\rho }\right)g\overline{\rho }.`$ (8)
This transformation redistributes the total energy between $`E_{\mathrm{pair}}`$ and $`E_{\mathrm{embed}}`$. When

$$g=\frac{dF}{d\overline{\rho }}|_{\overline{\rho }=\overline{\rho }_{\mathrm{xtal}}},$$
(9)
$`F^{}(\overline{\rho })`$ has a minimum at $`\overline{\rho }_{\mathrm{xtal}}`$, where $`\overline{\rho }_{\mathrm{xtal}}`$ is the value of $`\overline{\rho }`$ in the equilibrium crystal. This choice is called the effective pair format, and has been suggested as the most natural way to partition the energy between the pair and many-body contributions.Johnson and Oh (1989) In this format, when $`\overline{\rho }=\overline{\rho }_{\mathrm{xtal}}`$ the pair potential controls the energy change for any change of configuration that does not significantly alter $`\overline{\rho }`$, hence the name. More specifically, by performing a Taylor expansion about this reference density, one can show that to first order, the change in energy is due to the pair potential. Consequently, it is also a much more helpful format for relating the structure to the form of the potential.

The Gupta potential in this effective pair format becomes

$`\varphi _{\mathrm{eff}}\left(r\right)`$ $`=`$ $`2Ae^{p(r/r_01)}{\displaystyle \frac{\xi }{\sqrt{\overline{\rho }_{\mathrm{xtal}}}}}e^{2q(r/r_01)}`$ (10)
$`F_{\mathrm{eff}}\left(\overline{\rho }\right)`$ $`=`$ $`\xi \sqrt{\overline{\rho }}\left(1{\displaystyle \frac{1}{2}}\sqrt{{\displaystyle \frac{\overline{\rho }}{\overline{\rho }_{\mathrm{xtal}}}}}\right).`$ (11)
$`\varphi _{\mathrm{eff}}(r)`$ is a sum of two exponentials, so for $`p>2q`$ it is repulsive at short range and has an attractive well. The minimum in $`\varphi _{\mathrm{eff}}(r)`$ is at

$$r_{\mathrm{min}}=r_0\left(1+\frac{1}{p2q}\mathrm{log}\left(\frac{Ap\sqrt{\overline{\rho }_{\mathrm{xtal}}}}{\xi q}\right)\right)$$
(12)
and is of depth

$$\varphi _{\mathrm{eff}}\left(r_{\mathrm{min}}\right)=A\left(\frac{\xi q}{Ap\sqrt{\overline{\rho }_{\mathrm{xtal}}}}\right)^{p/p2q}\left(\frac{p2q}{q}\right)$$
(13)
As $`p2q`$ from above, the depth of the well in the pair potential goes to zero. $`F_{\mathrm{eff}}(\overline{\rho })`$ is quadratic in $`\sqrt{\overline{\rho }}`$ and has a minimum of depth $`\xi \sqrt{\overline{\rho }_{\mathrm{xtal}}}/2`$ at $`\overline{\rho }=\overline{\rho }_{\mathrm{xtal}}`$.

The Gupta parameters for the Zn and Cd potentials are given in Table 1, and the functions $`\varphi _{\mathrm{eff}}`$, $`\rho `$ and $`F_{\mathrm{eff}}`$ are shown in Fig. 1. The shallowness of the attractive well in the effective pair potential of cadmium is particularly apparent—it is only 0.26% of the depth of the minimum in the embedding function. The pair well depth is somewhat larger for zinc, but it is still only 0.70% of the embedding function minimum. For both systems the shallowness arises because $`p/q`$ is close to 2.

This feature of the potentials will have similar structural consequences for both systems. Firstly, the majority of the binding energy will come from the many-body interactions. Secondly, the pair potential provides relatively little constraint on the pair distances, except for the repulsion at short distances. Therefore, the most important feature for obtaining a low-energy configuration is to have the individual $`\overline{\rho }_i`$ values as close as possible to the optimal value. By contrast, there is little energetic advantage in having the nearest-neighbour distances close to the minimum of the pair potential. For a cluster the former can sometimes be more easily achieved through a disordered structure than one based on a lattice, because the additional flexibility of not having well-defined nearest-neighbour distances makes it easier to obtain close to optimal $`\overline{\rho }`$ values for the surface atoms. This is the source of the tendency to disorder for zinc and cadmium clusters described by the current potentials.

For bulk both zinc and cadmium are hexagonal close-packed with $`\overline{\rho }_{\mathrm{xtal}}^{\mathrm{Cd}}=8.042`$, and $`\overline{\rho }_{\mathrm{xtal}}^{\mathrm{Zn}}=8.638`$. These are anomalously low values for close-packed materials. For example, for an ideal close-packed crystal with all nearest-neighbours at $`r_0`$, the nearest-neighbour contribution to $`\overline{\rho }`$ is 12. However, the ratios of the unit cell parameters $`c/a`$ for the Zn and Cd crystals are particularly large, and so the contribution to $`\overline{\rho }`$ from nearest-neighbour distances with a component in the $`c`$ direction is significantly reduced.

The low values of $`\overline{\rho }_{\mathrm{xtal}}`$ have structural consequences for the clusters. An atom can achieve the optimal $`\overline{\rho }`$ values with only eight or nine neighbours, so this is possible for surface atoms. By contrast, the nearest-neighbour distances for atoms in the interior need to be elongated to prevent unfavourably large values of $`\overline{\rho }`$. This is bad news for cluster structures, such as the Mackay icosahedron, where the interior distances are naturally shorter than those on the surface.

The exponential nature of $`\rho (r)`$ means that significant changes to the $`\overline{\rho }`$ values can be achieved by relatively small changes to the nearest-neighbour distances. Consequently, it is somewhat easier for atoms, even those with low coordination numbers, to obtain nearly optimal $`\overline{\rho }`$ values. This behaviour contrasts with other many-body potential, for instance those produced by the force-matching method,Ercolessi and Adams (1994) where there is no presumed form for $`\rho (r)`$, and so more long-ranged functions can result.Lim et al. (1992); Ercolessi and Adams (1994); Mishin et al. (1999)

The current potentials were obtained by keeping the $`c/a`$ value fixed at the experimental value. This is at some cost in terms of the quality of the fit to other quantities. This led Cleri and Rosato to also construct potentials in which the $`c/a`$ ratio was allowed to vary in the fitting procedure, and resulted in a better quality fit for many properties.Cleri and Rosato (1993) Such a potential is available for cadmium and has dramatically different structural properties because $`p/q=3.49`$.

### II.2 Global Optimization

The global optimization of the zinc and cadmium clusters was performed using the basin-hopping Wales and Doye (1997); Wales and Scheraga (1999) (or Monte Carlo minimization Li and Scheraga (1987)) approach. This method has proved particularly successful in locating putative global minima for a wide variety of cluster systems. The nature of the potentials for the current clusters, particularly that $`p2q`$ is small, makes global optimization difficult compared to many other metal clusters of similar size. Therefore, a considerable computational effort was required to extend the results up to 125 atoms. It proved particularly important to supplement the standard unbiased runs from random starting points for each size, with short runs started from low-energy minima of nearby sizes with the appropriate number of atoms added or removed. As the structures of the two clusters are similar—seventy of the global minima are the same—it also proved useful to reoptimize the low-energy minima found for one metal for the other. These two types of ‘seeded’ runs were applied iteratively until no further new global minima were located.

It should be noted that, of course, there is no guarantee that I have been able to locate the true global minima, and the probability that a global minimum has been missed will increase with cluster size, as the size of the search space, and hence the number of minima,Tsai and Jordan (1993); Stillinger (1999); Doye and Wales (2002) increases exponentially with $`N`$. Examination of the statistics of how often independent runs locate the same lowest-energy minimum and the importance of the seeded runs can provide one with an idea of the likelihood that the true global minimum has been found. For example, for zinc clusters with less than ninety atoms virtually all the putative global minima were located in unbiased runs, but for $`N>100`$ the majority of putative global minima were only located in seeded runs. Similarly, for the cadmium clusters virtually all the putative global minima were only located in seeded runs for $`N65`$. This greater difficulty is because $`2pq`$ is closer to zero for cadmium. However, given the similarity of the observed structures for the two metals, I am confident that the vast majority of the putative global minima up to $`N=100`$ cannot be bettered, but beyond this size one’s degree of scepticism about the success of the global optimization should increase rapidly.

## III Global minima

The energies and point groups for the putative global minima are given in Tables 2 and 3. Point files are available online at the Cambridge Cluster Database.Web  The energies of the global minima are represented in Fig. 2 in such a way that makes particularly stable clusters stand out.

The previous results for these clusters were for a selection of sizes that often show highly symmetric structures, namely $`N=13`$, 38, 55 and 147.Michaelian et al. (2002) For $`N`$=13 and 38 the three lowest-energy minima reported in Ref. Michaelian et al., 2002 agree with the current results. However, for Zn55, Zn75 and Cd75 the lowest-energy structures reported by Michaelian et al. correspond at best to the third, eighth, and twenty-second lowest-energy minima, respectively, and lie 0.0031, 0.0186 and 0.0086 eV above the lowest-energy minima reported here. These energies are significant compared to the variations shown in Fig. 2. I did not systematically attempt to optimize 147-atom clusters, as locating the true global minimum for this size would be extremely difficult. However, short basin-hopping runs did find structures that were 0.0966 and 0.0304 eV lower in energy for Zn147 and Cd147, respectively, than the lowest-energy structures found by Michaelian et al.Michaelian et al. (2002)

Before, we examine the observed structures, there are a number of interesting features evident from Fig. 2. The energy zero in these figures is $`E_{\mathrm{ave}}`$, a four-parameter fit to the energies of the global minima, where the first two terms correspond to volume and surface energies. For these two clusters the surface term is exceptionally small. The ratio of the surface to the volume coefficient is 7.8% for Zn and 3.3% for Cd. For comparison the value of this ratio is 197% for Lennard-Jones clusters,Doye et al. (1999) 48% for Gupta lead clusters,Lai et al. (2002); Doy  60% for aluminium clusters,Doye (2003) and 93%, 87% and 63% for Sutton-Chen silver, nickel and gold clusters, respectivelyDoye and Wales (1998). One expects the surface energy to be lower for metal clusters than for a cluster interacting with a pair potential, because the lower coordinate surface atoms of a metal can increase their many-body embedding energy by shortening their nearest neighbour distances. But the Zn and Cd surface energies are very low even for metals, because the optimal value of $`\overline{\rho }`$ is small enough to be achievable by the surface atoms and the effective pair potential only makes a small contribution to the energy.

The second interesting feature is the magnitude of the fluctuations about the average energy. For Cd and Zn they are again exceptionally small. The average deviation from $`E_{\mathrm{ave}}`$ compared to the average energy per atom is 0.67% for Zn and 0.31% for Cd. For comparison, the value is 16% for Lennard-Jones clusters, 2.1% for Gupta lead clustersLai et al. (2002); Doy  and 5.6% for aluminium clustersDoye (2003) in the same size range.

The implications of these small fluctuations are that the differences between the more and less stable sizes indicated by Fig. 2 are small, making it more difficult to observe any such ‘magic’ numbers. Interestingly, one would expect that the properties of fully disordered clusters would evolve smoothly with size. By this measure, these zinc and cadmium clusters seem to be close to this limit.

We will first examine the zinc clusters in detail, and then later look at the relatively small structural difference between the two systems. A selection of zinc clusters are depicted in Fig. 3. that are either particularly stable or have some interesting structural feature. Up to $`N=10`$ the structure of the smallest zinc clusters are typical of what one usually finds for clusters modelled by empirical potentials. However, instead of structures leading up to the 13-atom icosahedra, more open structures are then preferred. The example shown in Fig. 3 for Zn13 can be considered as a polytetrahedral fragment of a 19-atom double icosahedron. It is insightful to examine, why this structure is lower in energy than the 13-atom icosahedron. The icosahedron is in fact $`0.335`$eV higher in energy, and this stems from an unfavourable embedding energy. The better pair energy, because of the icosahedron’s greater average coordination number, does little to offset this. As the distances between the central atom and the vertices of the icosahedron are 5% shorter than the distance between adjacent vertices and $`\rho (r)`$ increases rapidly with decreasing $`r`$, the central atom has an extremely high value of $`\overline{\rho }`$, namely 19.568, that is much greater than the optimal value. By contrast, the $`\overline{\rho }`$ values for all the atoms in the global minimum are relatively close to $`\rho _{\mathrm{xtal}}`$. This is achieved by the atoms with low coordination number having shorter average distances.

This example illustrates two important features of the potential. First, low-coordinate surface atoms are not disfavoured, in contrast to a potential with a strong pair component. Second, it is important that interior atoms do not have short distances. As we will see, it is often much better for an interior atom to have a coordination number larger than twelve since then its nearest-neighbour distances will be longer than those between its neighbours. This latter feature is illustrated by the 14-atom global minimum. The 13-coordinate atom has an average nearest-neighbour separation of $`1.052r_0`$, whereas the average separation for the nearest-neighbour contacts not involving this atom is $`0.958r_0`$. Consequently $`\overline{\rho }`$ for the central atom is only 10.381

The global minima for $`N`$=16–21 can all be characterized as distorted decahedral structures. Decahedral structures are based on pentagonal bipyramids (hence the name) and have a single five-fold axis of symmetry. The pentagonal bipyramids themselves are usually not particularly stable because they have a relatively non-spherical shape. By the introduction of reentrant grooves at the five equatorial vertices that are parallel to the five-fold axis, more stable Marks decahedra can be produced.Marks (1984); als  An 18-atom example is shown in Fig. 4 that was derived from a 23-atom pentagonal bipyramid, where the introduction of the grooves gives rise to five four-coordinate capping atoms.

This 18-atom Marks decahedron is in fact itself not that unstable and is only 0.042 eV above the global minimum. Even though the central atom is 12-coordinate, this atom is able to maintain a reasonable $`\overline{\rho }`$ value (9.867), whilst an expansion of the structure along the axial direction allows the surface distances to contract.

However, the energy can be further improved by distorting this Marks decahedron. The distortion of this decahedra is easy to see for the 17-atom global minimum. This structure can be derived from the Marks decahedron (minus one capping one atom) by a diamond-square-diamond rearrangement,Lipscomb (1966) where the contact between the atoms in the groove opposite the missing cap atom is broken and a contact between the adjacent capping atoms is formed (compare the second views of the two structures in Fig. 3 and 4). During this process the two capping atoms are also drawn closer to the central atom to form new contacts, thus giving the central atom a coordination number of 14. The planes containing the two five-fold rings of atoms that were parallel in the Marks decahedron are now splayed apart (see the third view of the Zn17 global minimum). It is also noticeable that the structure in the region of the distortion locally resembles the 13-atom icosahedron with the surface atoms now six coordinate.Soler et al. (2001) The second view of Zn17 looks down a two-fold axis of this icosahedral-like region.

The energetic advantage provided through this ‘groove-bridging’ distortion is mainly through the lowering of the energies of the two capping atoms that come into contact. The increase in their coordination number to 6 gives rise to both an increase in $`\overline{\rho }`$ and a more favourable pair energy. The energy of the central atom is left relatively unaffected, because the increase in coordination number is compensated by an increase in the average nearest-neighbour separation for that atom.

The other global minima for $`N`$=16–21 can be also understood in terms of this distortion. Zn16 simply has one less capping atom than Zn17. For $`N`$=18 and 19 two of the grooves are bridged, giving rise to a 16-coordinate central atom. The combined effect of the two distortions also causes the contact in the groove adjacent to the these two distortions to be broken, thus giving rise to a square face, which is itself capped in Zn19, and a four-fold axis of symmetry (see the second view of Zn19). Zn20 and Zn21, like Zn17, are only distorted at one groove; the additional atoms fill in some of the other grooves taking the structure closer towards the 23-atom pentagonal bipyramid.

Such distorted decahedral structures have previously been seen for a number of metal potentials,Doye and Wales (1998); Lloyd and Johnston (1998, 2000); Wilson and Johnston (2000); Soler et al. (2001); Darby et al. (2002); Lai et al. (2002); Sebetci and Güvenc (2003) particularly those that have a tendency to disorder. Indeed, Soler et al. have identified such structures as particularly important for amorphous gold clusters, and have used them as a basis to suggest schemes to produce potentially new magic number clusters for larger sizes.Soler et al. (2001)

As the size of the clusters increase, so of course must the number of atoms in the interior. Zn25 is typical of most of the clusters with two internal atoms. It can be considered as two interlocking distorted decahedra. In the view shown in Fig. 3 a Zn18-like fragment can be clearly seen. The exception to these type of structures occurs at Zn28. The six-fold symmetric structure is a continuation of Zn14 and has four stacked hexagonal rings arranged in an antiprismatic fashion. One would expect such a structure to have three internal atoms, but instead there is a vacancy in the very centre of the structure. Again, this is to avoid internal atoms with values of $`\overline{\rho }`$ that are too large.

Zn30 and Zn34 are typical of structures with three and four internal atoms. The overall shape reflects the triangular and tetrahedral arrangement of the internal atoms, but it is again hard to detect any overall order. Zn36 provides an interesting exception. The surface atoms have basically the same geometry as a recently identified 38-atom polytetrahedral structure.Doye et al. (2001, 2003) However, there are only four (not six) internal atoms, and these are arranged as a planar rhombus.

First at $`N`$=35 and then from $`N`$=38 a new series of distorted decahedral structures begins. Again they are based on a Marks decahedra that has been created from a pentagonal bipyramid (in this case with 54 atoms) by introducing reentrant grooves that are one layer deep. A similar mechanism of distortion is again seen. Atoms that are on the equatorial edges of the cluster either side of a groove are drawn closer together, although in this case a contact is not actually formed—the distortion is a diamond-square rather than a diamond-square-diamond process. The clusters at $`N`$=35, 39, 41, 45, 47 and 49 in Fig. 3 provide examples of structures in this sequence. In most of the examples, two such distortions are present at adjacent grooves, and so the whole edge between these grooves is brought closer to the centre of the cluster, forming new contacts with the closest interior atoms. As can be seen for the second view of the complete 49-atom Marks decahedron, this lead to a splaying of the planes of pentagonal rings that would otherwise have been parallel.

For these larger clusters, these distortions represent a smaller perturbation on the overall structure than for $`N`$=16-21. Indeed, they are barely visible when the cluster is viewed down the quasi-fivefold axis. All that can be seen is a slight opening of the angle at the reentrant groove.

Comparing the energetics of the distorted and undistorted 49-atom Marks decahedron shows that the main improvement arising from the distortion is for the interior atoms. The movement of the equatorial edges involved in the distortion closer to the centre of the cluster compresses the rest of the surface somewhat and increases the coordination number of some of the interior atoms. This leads to a structure where the ratio of the nearest neighbour distances for the interior atoms to those of the surface atoms increases from 1.049 to 1.060. Thus, there is less need for the surface to shrink inwards and compress the core in order to improve the $`\overline{\rho }`$ values of the surface atoms. Consequently, $`\overline{\rho }_{\mathrm{bulk}}`$, the average $`\overline{\rho }`$ value for the interior atoms, decreases from 9.197 to 8.756, i.e. closer to $`\overline{\rho }_{\mathrm{xtal}}`$, and so the energy of the interior atoms decreases. By contrast the distortion has little overall effect on the energy of the surface atoms. Generally, the atoms in the equatorial plane, especially those close to the distortion, improve their energies, but many of the other surface atoms lose out due to the breaking or stretching of contacts in the axial direction.

This decahedral series of structures is interrupted at N=42 and 43 by structures with four-fold symmetry. These clusters are loosely related to the 44-atom fcc octahedron, but with one or two opposite vertices removed. The top half of the cluster is then twisted with respect to the bottom to give a slightly buckled outer layer.

Beyond $`N`$=50 these groove-bridged Marks decahedra are no longer most stable. Instead, there is a size range where the structures typically have no overall order, but where motifs that resemble fragments of Marks decahedra and Mackay icosahedra are evident. Zn54 provides a typical example, whereas Zn52 is a structure with more order evident. From one side Zn52 has perfect Mackay icosahedral order. What has occurred to the other side can be determined by comparing to the 50-atom incomplete Mackay icosahedron shown in Fig. 4. The two contacts that complete the five-fold rings around the two empty vertices of the Mackay icosahedron have been broken, removing some of the tension in the outer layer of the Mackay icosahedron that makes it so energetically unfavourable.hiM  Then two atoms (oriented horizontally with respect to the views in the figures) are inserted into the coordination shell of the central atom, increasing its coordination number. Such a structure has previously been located for Au52 modelled by a Sutton-Chen potential.Doye and Wales (1998)

From the particularly stable structures based on the 18-atom and 49-atom Marks decahedra, it is not surprising that there is another series of magic numbers leading to the 100-atom Marks decahedron. Again this Marks decahedron is formed from a pentagonal bipyramid by the introduction of grooves of depth one layer. More generally, one can use this trend to predict potential magic number clusters outside the size range of this study. The sizes for such complete Marks decahedra are given by

$$N=\frac{5}{6}n^3+5n^2+\frac{61}{6}n+2$$
(14)
where $`n`$ is the number of atoms on an equatorial edge of the Marks decahedron. This gives $`N`$=18, 49, 100, 176, 282, $`\mathrm{}`$

The first of this set of groove-bridged decahedral structures occurs at $`N`$=75 and the last at $`N`$=101. Representative examples at $`N`$=75, 77, 79, 83, 86–88, 92, 96 and 100 are shown in Fig. 3. At the smallest sizes the decahedra are still quite asymmetric. From Fig. 2 one can see that the particularly stable sizes occur at $`N=1004m`$, where $`m`$=0–3. These structures can be formed from Zn100 by the sequential removal of 4-atoms from each groove to give structures that have $`m`$ grooves that are two layers deep. It is interesting to know that the complete Marks decahedron is not in fact the most stable of these structures, but that a slight asymmetry is preferred.

Comparing the energetics of the 100-atom distorted and undistorted Marks decahedra reveals a similar story to Zn49. The main stabilization of the distorted structure is due to the decrease in $`\overline{\rho }_{\mathrm{bulk}}`$ from 8.875 to 8.789. It is also interesting to analyse the reasons for the stability of structures at $`N`$=88 and 92 with four-coordinate surface atoms. Firstly, this is because a four-coordinate atom can compensate for its low coordination by having very short nearest-neighbour distances, thus achieving a reasonable $`\rho `$ value; e.g. 7.519 for Zn92. Secondly, this additional short contact reduces the need for the four surface atoms in contact with the adatom to shrink inward and compress the adjacent interior atoms. For Zn92 the $`\overline{\rho }`$ values for these atoms only increase by 0.313, despite this extra contact, because the pair separations for the other six contacts have expanded by 2.5%.

The structures at $`N`$=63-70 are also closely based on the 100-atom Marks decahedron, but with a further distortion. From the first view of Zn64 one can see the resemblance to the 63-atom structure in Fig. 4 that is an asymmetric fragment of the 100-atom distorted Marks decahedron. The second view of Zn64 shows the 2-fold axis. From the equivalent view of the 63-atom structure, it is clear that an extra atom has been added to the column of three capping atoms (at the bottom of the first view and in a vertical line in the centre of the top half of the cluster in the second view), and then a small twist has been given to the two halves of the structure. Zn69 has basically the same structure but with an increase in the size of one of the decahedral faces.

Based on the favourability of small distorted Marks decahedra (similarly to those we see at $`N`$=16–21) for gold clusters modelled by a Gupta potential, Soler et al. suggested a means for using this structural pattern to generate potentially stable large clusters.Soler et al. (2001) A 55-atom Mackay icosahedron can be generated from a 13-atom icosahedron by the addition of atoms at the centre of each nearest-neighbour contact. By the same process, an 80-atom structure can be generated from an 18-atom distorted Marks decahedron with one of the grooves bridged. Soler et al. envisaged that this structure would be most stable when the three four-coordinate atoms in this 80-atom structure were removed.Soler et al. (2001) However, the suggested structure is actually the global minimum for Zn80, reflecting the stabilization of low-coordinate atoms that is possible for this potential.

An alternative way to generate the Zn80 structure is to introduce reentrant grooves two layers deep into the 100-atom pentagonal bipyramid. This Marks decahedron is then distorted to produce four new contacts that bridge the two capping square pyramids adjacent to a groove. This distortion is analogous to that for bridging a groove one layer deep. In the region of the distortion, the structure looks locally like a 55-atom Mackay icosahedron.

Bridging a groove that is two layers deep is generally less favourable than bridging one that is one layer deep. There are only five examples for $`N`$100. The others occur for Zn71-73 and Zn84. Zn71 and Zn84 are quite interesting because they involve distortions of grooves that are both one and two layers deep. For example, Zn84 can be formed from Zn80, by first filling in one of the grooves, so that it is then only one layer deep and then applying the distortion to this groove.

After the completion of the previous decahedral series associated with the Marks decahedra at $`N`$=18 and 49, there are size ranges where the majority of clusters have no discernible overall order at $`N`$=22–34 and $`N`$=51–62, before structures based on the next Marks decahedron become lowest in energy. However, even though the next complete Marks decahedron is at $`N`$=176, beyond $`N`$=101 structures based on this larger Marks decahedron are immediately lowest in energy. From $`N`$=102–117 (excepting $`N`$=108) the zinc clusters are asymmetric groove-bridged decahedra with the same twist distortion as that for $`N`$=63-70. Like for Zn64 this leads to structures with a two-fold axis of symmetry at $`N`$=111 and 115 (see the second views of these structures in Fig. 4).

The first global minima that are fragments of the 176-atom Marks decahedron with only distortions bridging grooves one layer deep occur at $`N`$=122 and 123. There are also a number of structures, where grooves two layers deep are bridged; e.g. $`N`$=108, 120, 121 and 125. The structures at $`N`$=108, 121 and 125 are slightly more complicated, because the reentrant $`\{111\}`$ faces in the grooves of the undistorted Marks decahedra are not triangular but instead have two atoms along the outer apical edge. This is illustrated for Zn108, which is based on the 111-atom Marks decahedron shown in Fig. 4. On bridging this groove, it is favourable to remove three atoms that lie on the $`\sigma _h`$ mirror plane of the Marks decahedron, and to twist the top and bottom of the cluster slightly to remove a small $`\{100\}`$-type face that would otherwise result.

As has already been noted, for just over half the cadmium clusters in the size range considered, the global minimum have the same structures as the zinc clusters. A selection of examples where the global minimum differs from zinc are depicted in Fig. 5. The global minima are again dominated by distorted Marks decahedra, although the positions of capping atoms (e.g. Cd16 and Cd19) or the number of grooves that are bridged (e.g. Cd100) may be different. To take Cd100 as an example, this structure is second lowest in energy for zinc, because, although it has a better embedding energy than the global minimum, this is more than offset by a worse pair energy. However, because of the reduced magnitude of the pair interactions for cadmium this structure is now lowest in energy.

Another difference between the two systems is the more pronounced stability of the twisted Marks decahedra for cadmium, as indicated by Fig. 2 and the larger size ranges ($`N`$=63–72 and 102–121) for which these structures are most stable. A further difference is that decahedral structures where grooves of depth two are bridged are less common—there are only three examples for cadmium clusters in the size range considered here.

The least coincidence between the zinc and cadmium global minima probably occurs in the size windows between the series of distorted Marks decahedra where the structures often have no overall order. This probably reflects the large number of disordered structures that only have small differences in energy. Examples of clusters from these size ranges with $`N`$=26, 30, 38, 55 and 56 are shown in Fig. 5.

## IV Conclusion

I have analysed the structures of the global minima up to $`N`$=125 for two metallic potentials that have been previously found to have a particularly strong tendency to exhibit disordered clusters. This study confirms that the clusters exhibit none of usual ordered forms for materials that are close-packed in bulk. Instead, the majority of the clusters are based on distorted oblate Marks decahedra, but where the distortions are well-defined. There are series of structures associated with the 18-, 49-, 100- and 176-atom complete Marks decahedra. However, in between there are size windows where the majority of the clusters have no overall order, except for $`N>100`$ where there is a direct transition between structures based on the 100- and 176-atom Marks decahedra.

Because the effective pair potentials of these zinc and cadmium potentials are very shallow, thus providing relatively little constraint on the pair separations, there is a strong tendency for the surface atoms to contract inwards in order for these atoms to obtain a better many-body embedding energy. The resulting compression of the interior of the cluster causes the conventional structures to be disfavoured,Soler et al. (2000) because it results in a $`\overline{\rho }_{\mathrm{bulk}}`$ that is significantly larger than $`\overline{\rho }_{\mathrm{xtal}}`$. Instead structures for which the nearest-neighbour distances for the surface atoms are naturally longer than those for the interior atoms are likely to be favoured, because this allows the difference between $`\overline{\rho }_{\mathrm{bulk}}`$ and $`\overline{\rho }_{\mathrm{surf}}`$ to be reduced, bringing them both closer to the optimal value. The distortions of the oblate Marks decahedra help to achieve just this.

This paper should be seen in the context of a growing research program that has sought to understand how the observed structures for metal clusters depend on the form of the many-body potential.Soler et al. (2000); Baletto et al. (2002); Hendy and Doye (2002); Doye and Hendy (2003); Doye (2003) Even for simple potentials of the embedded-atom (or glue) form (Eq. 2) this leads to significant complexity. The current results further highlight how many-body potentials can lead to the stabilization of unusual structural forms.Soler et al. (2000, 2001); Hendy and Doye (2002); Doye and Hendy (2003); Doye (2003) Of course, at sufficiently large sizes the bulk structure, in this case hexagonal close-packed, should become most stable, but the energy differences found in this study suggest that the current clusters are far from this limit.

The emphasis in this paper has been on the zinc and cadmium clusters as model systems with a strong tendency to disorder. As noted in Sect. II.1, with the Gupta potential it is difficult to capture the large $`c/a`$ ratio for Zn and Cd without introducing discrepancies for other properties. I have also performed a brief set of global optimization runs using Cleri and Rosato’s alternative Gupta potential for Cd, where the $`c/a`$ ratio was allowed to vary during the fitting. For the selection of sizes I tested the usual ordered forms, i.e. close-packed, icosahedral and decahedral, were always most stable. Reyes-Nava et al. also found these cadmium clusters to have strong first-order-like melting transitions, again reflecting the strong ordering for this alternative potential.Reyes-Nava et al. (2003) Given this strong dependence on the potential parameterization, one should be somewhat sceptical about the abilities of these potentials to provide realistic models for zinc and cadmium clusters. However, the tendency to disorder for these clusters seems to be reproduced by electronic structure calculations using density-functional theory.Michaelian et al. (2002)

There has been little other relevant work. Experiments on these clusters have focussed on their electronic shell structure, Katakuse et al. (1986); Ruppel and Rademann (1992) rather than their geometries, and their are only a few theoretical studies to go beyond very small sizes. For cadmium clusters electronic structure calculations have been performed up to 20 atoms.Zhao (2001) These found Cd13 to be an icosahedron and the larger clusters to be based upon this structure.Zhao (2001) Ramprasad and Hoagland studied zinc clusters using a many-body potential of the embedded-atom form (Eq. 2), but where the effective interactions looked very different from the current Gupta potential.Ramprasad and Hoagland (1993) Of the series of candidate structures that they reoptimized, they also found icosahedral structures to be lowest in energy.

###### Acknowledgements.

The author is grateful to the Royal Society for the award of a University Research Fellowship.


format: text

text:
# Asymmetric Landau-Zener tunneling in a periodic potential

## Abstract

Using a simple model for nonlinear Landau-Zener tunneling between two energy bands of a Bose-Einstein condensate in a periodic potential, we find that the tunneling rates for the two directions of tunneling are not the same. Tunneling from the ground state to the excited state is enhanced by the nonlinearity, whereas in the opposite direction it is suppressed. These findings are confirmed by numerical simulations of the condensate dynamics. Measuring the tunneling rates for a condensate of rubidium atoms in an optical lattice, we have found experimental evidence for this asymmetry.

The phenomenon of Landau-Zener (LZ) tunneling landzen  is a basic quantum mechanical process. It is based on the solution of the Schrödinger equation for a two-level dynamics when a parameter of the Hamiltonian system is time dependent. If at time $`t=\mathrm{}`$ the system is prepared in one adiabatic state of the Hamiltonian, the time dependence of the Hamiltonian implies that at time $`t=+\mathrm{}`$ there is a finite probability that the system will occupy the other adiabatic state. As far as this tunneling behaviour is concerned, complete symmetry exists between the adiabatic states. Variations of the LZ model have been studied mullen-ben-tzt ; gefen-ben-zt ; fish-mul-benj , and an observation of LZ dynamics in classical optical systems has been reported bou-dekker-lzdyn . More recently, LZ tunneling within a periodic potential was studied for a nonlinear two-level system in which the level energies depend on the occupation of the levels zobay ; wuniu-lz . It was discovered that a nonlinearity with a positive sign enhances the tunneling probability between the ground band and the first excited band. Moreover, Niu and coworkers discovered a nonzero LZ tunneling probability even in the fully adiabatic limit when the nonlinearity was larger than a critical value wuniu-lz . Critical values for deformations of the energy level structures were obtained in refs. zobay ; pethick ; mueller . In a Bose-Einstein condensate inside a periodic potential such as an optical lattice, the mean-field interaction between the atoms can be comparable to other energy scales of the system and hence the level-dependent energy shift can lead to an observable modification of the tunneling behaviour.

In the present work we explore, theoretically and experimentally, the Landau-Zener tunneling between Bloch bands of a Bose-Einstein condensate in an accelerated optical lattice. The optical lattice depth controls the tunneling barrier, while the optical lattice acceleration controls the time dependence of the Hamiltonian. We show that the mean-field nonlinearity produces an asymmetry for the tunneling probability between the adiabatic states of the Hamiltonian. More precisely the tunneling probability from the lower energy adiabatic state to the upper one is enhanced, while the tunneling probability from the higher energy level to a lower level is suppressed. Numerical integration of the one-dimensional Gross-Pitaevskii equation and a simple two-state model demonstrate this asymmetry. Moreover, our experimental data obtained with a rubidium Bose-Einstein condensate confirms this prediction.

The qualitative explanation of this asymmetry in the tunneling transition probabilities depends on the choice of the quantum representation used for the description of the nonlinear system. We have identified two different basic mechanisms, acting simultaneously or independently. In the first mechanism the nonlinear term of the Schrödinger equation acts as a perturbation whose strength is proportional to the energy level occupation. If the initial state of the condensate in the lattice corresponds to a filled lower level of the state model, then the lower level is shifted upward in energy while the upper level is left unaffected. This reduces the energy gap between the lower and upper level and enhances the tunneling. On the contrary, if all atoms fill the upper level then the energy of the upper level is increased while the lower level remains unaffected. This enhances the energy gap and reduces the tunneling. The overall balance leads to an asymmetry between the two tunneling processes. In a different representation the nonlinearity term in the Schrödinger equation produces an additional term to the optical lattice depth. This additional term increases or decreases the optical lattice depth depending on whether the tunneling process proceeds from the lower level to the upper one, or vice versa. Thus the tunneling due to a weaker or stronger potential barrier is increased or decreased, as the case may be, depending on the initial state.
The motion of a Bose-Einstein condensate in an accelerated 1D optical lattice is described by the Gross-Pitaevskii equation

$$\begin{array}{c}i\mathrm{}\frac{\psi }{t}=\frac{1}{2M}\left(i\mathrm{}\frac{}{x}Ma_Lt\right)^2\psi +\hfill \\ \hfill +\frac{V_0}{2}\mathrm{cos}(2k_lx)\psi +\frac{4\pi \mathrm{}^2a_s}{M}\left|\psi \right|^2\psi \end{array}$$
(1)
where $`M`$ is the atomic mass, $`k_L=\pi /d`$ is the optical lattice wavenumber with $`d`$ the optical lattice step, and $`V_0`$ is the strength of the periodic potential depth. The $`s`$-wave scattering length $`a_s`$ determines the nonlinearity of the system. Equation 1 is written in the comoving frame of the lattice, so the inertial force $`Ma_L`$ appears as a momentum modification. The wavefunction $`\psi `$ is normalized to the total number of atoms in the condensate and we define $`n_0`$ as the average uniform atomic density. Defining the dimensionless quantities $`E_{rec}=\mathrm{}^2k_L^2/2M`$, $`\stackrel{~}{x}=2k_Lx`$, $`\stackrel{~}{t}=8E_{rec}t/\mathrm{}`$, and rewriting $`\stackrel{~}{\psi }=\psi /\sqrt{n_0}`$, $`\stackrel{~}{v}=V_0/16E_{rec}`$, $`\stackrel{~}{\alpha }=Ma_L/16E_{rec}k_L`$, $`C=\pi a_sn_0/k_L^2`$, Eq. (1) is cast in the following form wuniu-lz :

$$\begin{array}{c}i\frac{\psi }{t}=\frac{1}{2}\left(i\frac{}{x}\alpha t\right)^2\psi +v\mathrm{cos}(x)\psi +C\left|\psi \right|^2\psi \hfill \end{array}$$
(2)
where we have replaced $`\stackrel{~}{x}`$ with $`x`$, etc. In the neighborhood of the Brillouin zone edge we can approximate the wave function by a superposition of two plane waves (the two level model of ref. wuniu-lz ), assuming that only the ground state and the first excited state are populated. We then substitute $`\psi (x,t)=a(t)e^{iqx}+b(t)e^{i(q1)x}`$, with $`|a(t)|^2+|b(t)|^2=1`$ in Eq. (2). Comparing the coefficients of $`e^{iqx}`$ and $`e^{i(q1)x}`$, linearizing the kinetic terms and dropping the irrelevant constant energy $`1/8+C[1+(|a|^2+|b|^2)/2]`$, Eq. (2) assumes the form

$$\begin{array}{c}i\frac{}{t}\left(\begin{array}{c}a\\ b\end{array}\right)=\left[\frac{\alpha t}{2}\sigma _3+\frac{v}{2}\sigma _1\right]\left(\begin{array}{c}a\\ b\end{array}\right)+\frac{C}{2}(|b|^2|a|^2)\sigma _3\left(\begin{array}{c}a\\ b\end{array}\right)\hfill \end{array}$$
(3)
where $`\sigma _ii=1,2,3`$ are the Pauli matrices note1 . The adiabatic energies of Eq. 3 have a butterfly structure at the band edge of the Brillouin zone for $`Cv`$ wuniu-lz ; pethick ; mueller , but in the present work we always work in a regime where $`Cv`$, hence that structure plays no role.

In the linear regime $`(C=0)`$, evaluating the transition probability in the adiabatic approximation, we find the linear LZ formula for the tunneling probability $`r`$

$$r=e^{\frac{\pi v^2}{2\alpha }}$$
(4)
expressing the occupation changes in terms of the rate $`\alpha `$ at which the diagonal energies of the linear Hamiltonian change their value, and of the off-diagonal interaction energy $`v`$. In the nonlinear regime, as the nonlinear parameter $`C`$ grows, the lower to upper tunneling probability grows as well until an adiabaticity breakdown occurs at $`C=v`$ wuniu-lz . The upper to lower tunneling probability, on the other hand, decreases with increasing nonlinearity note2 . We derived the tunneling rate from the numerical integration of Eq. (3). In Fig. 1(a) we plot the lower to upper tunneling rate (initial $`(a,b)`$ =$`(1,0)`$) and the upper to lower tunneling rate (initial $`(a,b)`$ =$`(0,1)`$) of the Bose-Einstein condensate as a function of the nonlinear parameter $`C`$ for different accelerations of the optical lattice. We see that for $`C=0`$ the rate is the same for both tunneling directions whereas for $`C0`$ the two rates are different, and the smaller the acceleration the larger the difference. This result is intuitive since for very small accelerations the main contribution originates from the nonlinear effect (the linear tunneling of Eq. (4) being small) while for large accelerations the main contribution comes from the linear effect. We confirmed the presence of a tunneling asymmetry by integrating directly Eq. (1) (taking into account the full experimental protocol described below), finding qualitative agreement with the prediction of the two-state model. For small $`C`$ values, we have fitted the $`C`$ dependence of the tunneling rate of Fig.1 (a) through the following expression zobay ; wuniu-lz :

$$r(C)=e^{\frac{\pi v^2}{2\alpha }(1\pm \beta \frac{C}{v})}$$
(5)
with $`\beta =0.75,0.17,0.14`$ for the different acceleration values.

The nonlinear regime is more easily interpreted by writing Eq. ((3)) as

$$\begin{array}{c}i\frac{}{t}\left(\begin{array}{c}a\\ b\end{array}\right)=\left[\frac{\alpha t}{2}\sigma _3+\frac{v}{2}\sigma _1\right]\left(\begin{array}{c}a\\ b\end{array}\right)\frac{C}{2}\left(\begin{array}{cc}|a|^2& b^{}a\\ a^{}b& |b|^2\end{array}\right)\left(\begin{array}{c}a\\ b\end{array}\right)\hfill \end{array}$$
(6)
Here the nonlinear diagonal terms represent a shift of the energy levels, different for an occupation of the ground or upper level, decreasing or increasing the energy gap depending on the initial state. The off-diagonal terms of Eq. 6 produce an equivalent contribution to the nonlinear tunneling process. The off-diagonal terms modify the interaction term $`v`$ equivalent to a Rabi frequency in the two-level model. The off-diagonal scalar product between the two states $`a^{}b`$, calculated applying the adiabatic approximation technique of crisp , changes sign depending on the initial state (lower or upper). Thus the linear Rabi frequency is modified by the nonlinear off-diagonal term, and the increase or decrease depends on the initial state. This nonlinear modification of the interaction energy constitutes an additional contribution to the asymmetric tunneling rate.

Experimentally, we investigated the phenomenon of asymmetric tunneling between the energy bands of Bose-Einstein condensates in an optical lattice using a setup described in detail in morsch01 ; cristiani02 . Briefly, we create condensates of $`N10^4`$ rubidium atoms in a time-orbiting potential (TOP) trap. Once condensation has been achieved, the mean trapping frequency $`\overline{\nu }_{trap}`$ of the magnetic trap is adiabatically reduced to values between $`15\mathrm{Hz}`$ and $`50\mathrm{Hz}`$. Thereafter, two laser beams with waists of $`1.8\mathrm{mm}`$ and intersecting at an angle $`\theta =38\mathrm{deg}`$ at the position of the condensate are switched on with a linear ramp of duration $`\tau _{ramp}=10\mathrm{ms}`$, thus ensuring adiabaticity of the loading process. The beams are detuned to the red side of the rubidium atomic resonance by $`30\mathrm{GHz}`$ and have a variable frequency difference $`\mathrm{\Delta }\nu `$ between them, controllable through two acousto-optic modulators which are also used to vary the intensity of the beams. In this way, a periodic potential with lattice constant $`d=1.18\mu \mathrm{m}`$ and lattice recoil energy $`E_{rec}/h=455\mathrm{Hz}`$ is created, which through the frequency difference $`\mathrm{\Delta }\nu `$ can be made to move at a constant velocity $`v=d\mathrm{\Delta }\nu `$ or accelerated with $`a_L=d\frac{d\mathrm{\Delta }\nu }{dt}`$. For the current experiment, lattice depths between $`0.25E_{rec}`$ and $`2.5E_{rec}`$ were used.

Landau-Zener tunneling between the two lowest energy bands of a condensate inside an optical lattice is investigated in the following way. Initially, the condensate is loaded adiabatically into one of the two bands. Subsequently, the lattice is accelerated in such a way that the condensate crosses the edge of the Brillouin zone once, resulting in a finite probability for tunneling into the other band (higher-lying bands can be safely neglected as their energy separation at the edge of the Brillouin zone is much larger than the band gap). After the tunneling event, the two bands have populations reflecting the Landau-Zener tunneling rate (assuming that, initially, the condensate populated one band exclusively). In order to experimentally determine the number of atoms in the two bands, we then increase the lattice depth (from $`2E_{rec}`$ to $`4E_{rec}`$) and decrease the acceleration (from $`3\mathrm{m}\mathrm{s}^2`$ to $`2\mathrm{m}\mathrm{s}^2`$). In this way, successive crossings of the band edge will result in a much reduced Landau-Zener tunneling probability between the ground state band and the first excited band (of order a few percent). The fraction of the condensate that after the first tunneling event populated the ground state band will, therefore, remain in that band, whereas the population of the first excited band will undergo tunneling to the second excited band with a large probability (around $`90`$ percent) as the gap between these two bands is smaller than the gap between the two lowest bands by a factor $`5`$ for our parameters. Once the atoms have tunneled into the second excited band, they essentially behave as free particles since higher-lying band-gaps are smaller still, meaning that the fraction of the condensate that populated the first excited band after the first tunneling event will no longer be accelerated by the lattice. In summary, using this experimental sequence we selectively accelerate that part of the condensate further that populates the ground state band. In practice, in order to get a good separation between the two condensate parts after a time-of-flight, we accelerate the lattice to a final velocity of $`46v_{rec}`$ footnote\_finvel  and absorptively image the condensate after $`22\mathrm{ms}`$ (see Fig. 2).

In order to investigate tunneling from the ground state band to the first excited band, we adiabatically ramped up the lattice depth with the lattice at rest and then started the acceleration sequence. The tunneling from the first excited to the ground-state band is investigated in a similar way, except that in this case we initially prepare the condensate in the first excited band by moving the lattice with a velocity of $`1.5v_{rec}`$ (through the frequency difference $`\mathrm{\Delta }\nu `$ between the acousto-optic modulators) when switching it on. In this way, in order to conserve energy and momentum the condensate must populate the first excited band at a quasi-momentum half-way between zero and the edge of the first Brillouin zone nist . Thereafter, the same acceleration sequence as described above is used in order first to induce Landau-Zener tunneling and then to separate the fraction of the condensate that underwent tunneling to the lowest band from the one that remained in the first excited band. For both tunneling directions, the tunneling rate is measured as

$$r=\frac{N_{tunnel}}{N_{tot}},$$
(7)
where $`N_{tot}`$ is the total number of atoms measured from the absorption picture. For the tunneling from the first excited band to the ground band, $`N_{tunnel}`$ is the number of atoms accelerated by the lattice, i.e. those detected in the final velocity class $`4v_{rec}`$, whereas for the inverse tunneling direction, $`N_{tunnel}`$ is the number of atoms detected in the $`v=0`$ velocity class.

To illustrate our method for measuring the tunneling in both directions, in Fig. 3 we report the two tunneling rates as a function of lattice depth for a condensate in a weak magnetic trap and hence a small value of the interaction parameter $`C`$ note3 . In this case, both tunneling rates are essentially the same and agree well with the linear Landau-Zener prediction. By contrast, when $`C`$ is increased, the two tunneling rates begin to differ, as can be seen in Fig. 1(b). Qualitatively we find agreement with the theoretical predictions of the non-linear Landau-Zener model, whereas quantitatively there are significant deviations. We believe these to be partly due to experimental imperfections. In particular, the sloshing (dipolar oscillations) of the condensate inside the magnetic trap can lead to the condensate not being prepared purely in one band due to non-adiabatic mixing of the bands if the initial quasimomentum is too close to a band-gap. Furthermore, a numerical simulation of the experiment shows that for large values of $`C`$, for which the magnetic trap frequency was large, the measured tunneling rates are significantly modified by the presence of the trap. We have, however, verified that when $`C`$ in the simulation is varied without varying the trap frequency, the asymmetric tunneling effect persists.
In summary, we have numerically simulated Landau-Zener tunneling between two energy bands in a periodic potential and found that, in the presence of a nonlinear interaction term, an asymmetry in the tunneling rates arises. Experimentally, we have measured these tunneling rates for different values of the interaction parameter and found qualitative agreement with the simulations. We note here that the phenomenon of asymmetric tunneling should be a rather general feature of quantum systems exhibiting a nonlinearity. For instance, calculating the energy shift due to a nonlinearity for two adjacent levels of a harmonic oscillator, one finds that both levels are shifted upwards in energy, the shift being proportional to the population of the respective level. The energy difference between the levels, therefore, decreases if only the lower state is populated and increases if all the population is in the upper level. Finally, we note that state-dependent mean-field shifts have also been observed in measurements of the clock-shift in ultra-cold and Bose-condensed atomic samples harber02 .
This work was supported by the MURST (COFIN2000 Initiative), the INFM (PRA ’Photonmatter’), and by the EU through Contract No. HPRN-CT-2000-00125.


text:
# Bohmian Quantum Gravity in the Linear Field Approximation

## I INTRODUCTION

There are several approaches to the quantization of gravity. Some important ones are, the canonical formalism (or the WDW formalism), Narlikar and Padmanabhan approach to quantization of the conformal degree of freedom of the space–time metric, the standard Bohmian quantum gravity, and a new approach relating quantum mechanics and the geometry of the space–time. Each of the above theories have some positive and some negative aspects, so that there is no universally accepted theory of quantum gravity. As some authors have mentioned, Bohmian quantum gravity has some noticable aspects. It has not the problem of time, the measurement problem and the difficulty with the meaning of the wavefunction. Also it presents a trajectory for the system just like in the classical domain.

To understand how this theory works, let us briefly explain the Bohmian quantum mechanics of a single non–relativistic particle. In the standard quantum mechanics, the particle’s wavefunction satisfies the Schrödinger equation:

$$i\mathrm{}\frac{\mathrm{\Psi }}{t}=\frac{\mathrm{}^2}{2m}^2\mathrm{\Psi }+V(\stackrel{}{x};t)\mathrm{\Psi }$$
(1)
If one decomposes the wavefunction into its norm and its phase:

$$\mathrm{\Psi }=R\mathrm{exp}(iS/\mathrm{})$$
(2)
one gets the following two equations:

$$\frac{S}{t}+\frac{|\stackrel{}{}S|^2}{2m}+V+𝓠=0$$
(3)
$$\frac{R^2}{t}+\stackrel{}{}\left(R^2\frac{\stackrel{}{}S}{m}\right)=0$$
(4)
where

$$𝓠=\frac{\mathrm{}^2}{2m}\frac{^2R}{R}$$
(5)
As it can be seen, we have a modified Hamilton–Jacobi equation and the continuity equation provided we assume that the particle’s trajectory is given by the classical relation:

$$\stackrel{}{v}=\frac{\stackrel{}{}S}{m}$$
(6)
The Newton’s equation of motion can be derived by taking the gradiant of the Hamilton–Jacobi equation and using the guidance relation (6). The result is:

$$m\frac{d^2\stackrel{}{x}}{dt^2}=\stackrel{}{}V\stackrel{}{}𝓠$$
(7)
So the Bohmian interpretation is as follows. To quantize any classical system add a quantum potential to the classical Hamilton–Jacobi equation. The quantum potential is given in terms of the density of an hypothetical ensemble of the system ($`R^2`$) and thus add the continuity equation to have a self–consistent system of equations.

Bohmian quantum gravity is also constructed successfully. But it has also some problems. The most essential one is perhaps the lack of general covariance. Although some authors look at this property as an essential character of quantum gravity.

Investigating Bohmian quantum gravity in the linear field approximation of gravity is perhaps a good idea, because as it is the case for classical gravity some important effects can be understood better, in the linear approximation. In this paper we shall do this and then study the theory for the Newtonian limit and the static spherically symmetric metric, and some observable phenomena are predicted.

## II LINEAR FIELD BOHMIAN QUANTUM GRAVITY

Our plan in this section is to make Bohmian quantum gravity in the linear field approximation, i.e. when the space–time metric is expanded about the flat Minkowski space–time as:

$$g_{\mu \nu }=\eta _{\mu \nu }+h_{\mu \nu }$$
(8)
and considering only linear terms in $`h_{\mu \nu }`$. The lagrangian density of the linear field classical gravity is given by:

$$𝓛=\frac{1}{2}_\lambda \varphi _{\mu \nu }^\lambda \varphi ^{\mu \nu }$$
(9)
where as it is convenient we introduced the new field variables:

$$\varphi _{\mu \nu }=h_{\mu \nu }\frac{1}{2}h\eta _{\mu \nu }$$
(10)
The canonical momenta are defined as:

$$𝚷_{\mu \nu }=\frac{\delta 𝓛}{\delta \dot{\varphi }^{\mu \nu }}=\frac{\varphi _{\mu \nu }}{t}$$
(11)
and the hamiltonian density is:

$$𝓗=\frac{1}{2}𝚷_{\mu \nu }𝚷^{\mu \nu }+\frac{1}{2}_i\varphi _{\mu \nu }^i\varphi ^{\mu \nu }$$
(12)
The general coordinate transformation invariance of the linear field gravity is guaranteed by choosing some specific gauge. The common gauge is the harmonic gauge:

$$^\mu \varphi _{\mu \nu }=0$$
(13)
The transition to the quantum theory can be achieved via Dirac’s canonical quantization scheme, in which we replace:

$$𝚷_{\mu \nu }i\mathrm{}\frac{\delta }{\delta \varphi _{\mu \nu }}$$
(14)
leading to the following wave equation:

$$i\mathrm{}\frac{\mathrm{\Psi }}{t}=\left(\frac{\mathrm{}^2}{2}d^3x^{}\frac{\delta ^2}{\delta \varphi _{\mu \nu }(x^{})\delta \varphi ^{\mu \nu }(x^{})}+\frac{1}{2}d^3x^{}_i^{}\varphi _{\mu \nu }^i\varphi ^{\mu \nu }\right)\mathrm{\Psi }$$
(15)
where $`\mathrm{\Psi }[\varphi _{\alpha \beta }]`$ is the wavefunctional.

The Bohmian interpretation of this wave equation can be achieved via setting:

$$\mathrm{\Psi }=R\mathrm{exp}(iS/\mathrm{})$$
(16)
leading to the Hamilton–Jacobi equation:

$$\frac{S}{t}+\frac{1}{2}d^3x^{}\frac{\delta S}{\delta \varphi _{\mu \nu }(x^{})}\frac{\delta S}{\delta \varphi ^{\mu \nu }(x^{})}+\frac{1}{2}d^3x^{}_i^{}\varphi _{\mu \nu }(x^{})^i\varphi ^{\mu \nu }(x^{})+𝓠=0$$
(17)
and the continuity equation:

$$\frac{R^2}{t}+d^3x^{}\frac{\delta }{\delta \varphi _{\mu \nu }(x^{})}\left(R^2\frac{\delta S}{\delta \varphi ^{\mu \nu }(x^{})}\right)=0$$
(18)
with the quantum potential defined as:

$$𝓠=\frac{\mathrm{}^2}{2R}d^3x^{}\frac{\delta ^2R}{\delta \varphi _{\mu \nu }(x^{})\delta \varphi ^{\mu \nu }(x^{})}$$
(19)
The Bohmian trajectories can be obtained via the guidance equation:

$$\frac{\varphi _{\mu \nu }}{t}=\frac{\delta S}{\delta \varphi ^{\mu \nu }}$$
(20)
and the gauge condition is given by the equation (13).

It is important here to note that we have not included the gauge condition as it is always done. That is to say we have not set:

$$_\mu \varphi ^{\mu \nu }=0\left(\mathrm{\Pi }^{0\nu }+_i\varphi ^{i\nu }\right)|\mathrm{\Psi }=\left(i\mathrm{}\frac{\delta }{\delta \varphi _{0\nu }}+_i\varphi ^{i\nu }\right)|\mathrm{\Psi }=0$$
(21)
Because this forces the state to satisfy the gauge condition and the Bohmian trajectories are not necessarily consistent with it. If one wishes to save the gauge condition at the level of trajectories, one should apply it on the Bohmian trajectories $`\varphi ^{\mu \nu }(\stackrel{}{x},t)`$.

The field equation can be derived by taking variation of the Hamilton–Jacobi equation (17) with respect to $`\varphi _{\mu \nu }`$. The result is:

$$\mathrm{}\varphi _{\mu \nu }(x)=\frac{\delta 𝓠}{\delta \varphi ^{\mu \nu }(x)}$$
(22)
Let us now, apply these results to the Newtonian approximation. In this limit the metric is given by:

$$h_{00}=h=2\mathrm{\Phi };h_{ij}=h_{0i}=0$$
(23)
where $`\mathrm{\Phi }`$ is the Newton’s gravitational potential. The $`\varphi `$–field is given by:

$$\varphi _{00}=\mathrm{\Phi };\varphi _{ij}=\mathrm{\Phi }\delta _{ij};\varphi _{0i}=0$$
(24)
and we have the identities:

$$\frac{\delta }{\delta \varphi _{00}}=\frac{\delta }{\delta \mathrm{\Phi }};\frac{\delta }{\delta \varphi _{ij}}=\delta _{ij}\frac{\delta }{\delta \mathrm{\Phi }}$$
(25)
So the Hamilton–Jacobi equation and the continuity equation are given by:

$$\frac{S}{t}+d^3x^{}\left(\frac{\delta S}{\delta \mathrm{\Phi }(x^{})}\right)^2+\frac{1}{2}d^3x^{}_i^{}\mathrm{\Phi }(x^{})^i\mathrm{\Phi }(x^{})+𝓠=0$$
(26)
$$\frac{R^2}{t}+4d^3x^{}\frac{\delta }{\delta \mathrm{\Phi }(x^{})}\left(R^2\frac{\delta S}{\delta \mathrm{\Phi }(x^{})}\right)=0$$
(27)
where the quantum potential is given by:

$$𝓠=\frac{2\mathrm{}^2}{R}d^3x^{}\frac{\delta ^2R}{\delta \mathrm{\Phi }(x^{})^2}$$
(28)
The guidance relation is:

$$\frac{\mathrm{\Phi }}{t}=\frac{\delta S}{\delta \mathrm{\Phi }}$$
(29)
and the condition of the Newtonian approximation is given by:

$$\frac{\mathrm{\Phi }}{t}=0$$
(30)
instead of the relation (13). As a result of this last equation we have:

$$\frac{\mathrm{\Phi }}{t}=\frac{\delta S}{\delta \mathrm{\Phi }}=0$$
(31)
so we have encountered with a pure quantum case***That is a case in which quantum potential balances the classical potential. So the quantum potential is not ignorable and we have no classical limit. and so we have:

$$E+\frac{1}{2}d^3x^{}_i^{}\mathrm{\Phi }(x^{})^i\mathrm{\Phi }(x^{})=𝓠$$
(32)
where we have assumed $`S/t=E`$, because in the Newtonian limit physical quantities should not depend on time. The continuity equation is satisfied identically.

The field equation can be derived by differentiating the Hamilton–Jacobi equation with respect to $`\mathrm{\Phi }`$. The result is:

$$^2\mathrm{\Phi }=\frac{\delta 𝓠}{\delta \mathrm{\Phi }(x)}$$
(33)
Until now we have encountered with the free linear field Bohmian quantum gravity. That is we have not considered the matter. If one wants to add the matter field semi–classicallyIn order to take into account the quantum effects of matter, one should start from the linear field lagrangian density with matter, and then proceed as in the previous section. The only difference would be that the matter quantum potential is also present in the field equations. We did not do so here, because in this work, the quantum effects of gravity are of consider., one should add merely the source term $`\kappa \rho `$ to the above equation, where $`\rho `$ is the matter density and $`\kappa `$ is the gravitational constant. So we have:

$$^2\mathrm{\Phi }=\kappa \rho +\frac{\delta 𝓠}{\delta \mathrm{\Phi }(x)}$$
(34)
## III STATIC SPHERICALLY SYMMETRIC SOLUTION

Let’s now look for static spherically symmetric solution of the equations of the previous section. We do this in the case of Newtonian and weak field approximations.

* First we work in the Newtonian approximation. Let’s choose $`R`$ as a packet around the classical solution of the gravitational field equation:

$$R=\mathrm{exp}\left\{\frac{\alpha }{2}\left[\mathrm{\Phi }(x)\mathrm{\Phi }^{(c)}(x)\right]^2\right\}$$
(35)
where $`\mathrm{\Phi }^{(c)}`$ is the classical field, i.e. the field satisfying the equation (34) without $`𝓠`$. This leads to the following quantum potential:

$$𝓠=\alpha \mathrm{}^2\left[1\alpha (\mathrm{\Phi }\mathrm{\Phi }^{(c)})^2\right]$$
(36)
so the field equation is given by:

$$^2\mathrm{\Phi }=\kappa \rho \alpha ^2\mathrm{}^2(\mathrm{\Phi }\mathrm{\Phi }^{(c)})$$
(37)
Let us look for the gravitational potential, when we have a point–like source:

$$\rho =M\delta (\stackrel{}{x})$$
(38)
We solve the equation (37) by setting:

$$\mathrm{\Phi }=\mathrm{\Phi }^{(c)}+f=\frac{\kappa M}{r}+f$$
(39)
so we have:

$$^2f=\alpha ^2\mathrm{}^2f$$
(40)
The solution is given by:

$$\mathrm{\Phi }=\frac{\kappa M}{r}+\underset{lm}{}\left[a_{lm}j_l(\alpha \mathrm{}r)+b_{lm}n_l(\alpha \mathrm{}r)\right]Y_{lm}(\theta ,\phi )$$
(41)
where $`a_{lm}`$ and $`b_{lm}`$ are constants.

Some points about this solution is of interest:

+ Because of the asymptotic behaviour of the spherical Bessel functions:

    $$j_l(x)\frac{1}{x}\mathrm{sin}\left(x\frac{l\pi }{2}\right)\text{as}x\mathrm{}$$
    (42)
    $$n_l(x)\frac{1}{x}\mathrm{cos}\left(x\frac{l\pi }{2}\right)\text{as}x\mathrm{}$$
    (43)
    one can see that this quantum solution also behaves as $`\frac{1}{r}`$ at large distances from the source.
    + It must be noted that the weak field approximation is only applicable for large $`r`$’s. So this solution does not tell anything about the event horizon or about the singularity.
    + If one sets $`M=0`$, so that we have no source at all, we have a pure quantum solution in which we have source–less gravity. This means that the quantum fluctuations of gravity can produce observable gravity.
    + Because of the oscillatory nature of the spherical Bessel functions, the gravitational potential is also oscillatory. So we have a set of periodic stable points, where the potential is minimum, and a test particle can be at rest at that locations (This is true also for the case $`M=0`$). It must be noted that, if we choose a sharp wave packet ($`\alpha \mathrm{}`$ in equation (35)), the frequency of such oscillations is very large.
* Let’s now investigate the static spherically symmetric metric in the linear field approximation. The relations (17) and (18) would be now:

$$\frac{1}{2}d^3x^{}\frac{\delta S}{\delta \varphi _{\mu \nu }(x^{})}\frac{\delta S}{\delta \varphi ^{\mu \nu }(x^{})}+\frac{1}{2}d^3x^{}_i^{}\varphi _{\mu \nu }(x^{})^i\varphi ^{\mu \nu }(x^{})+𝓠=0$$
(44)
and

$$d^3x^{}\frac{\delta }{\delta \varphi _{\mu \nu }(x^{})}\left(R^2\frac{\delta S}{\delta \varphi ^{\mu \nu }(x^{})}\right)=0$$
(45)
The Bohmian trajectories can be obtained via the guidance equation:

$$0=\frac{\varphi _{\mu \nu }}{t}=\frac{\delta S}{\delta \varphi _{\mu \nu }}$$
(46)
and the gauge condition (13). It is a simple task to see that the corrected field equation is given by:

$$^2\varphi _{\mu \nu }=\frac{\delta 𝓠}{\delta \varphi ^{\mu \nu }(x)}$$
(47)
As in the Newtonian case, choosing:

$$R=\mathrm{exp}\left\{\frac{\alpha }{2}\left[\varphi _{\mu \nu }\varphi _{\mu \nu }^{(c)}\right]\left[\varphi ^{\mu \nu }\varphi ^{(c)\mu \nu }\right]\right\}$$
(48)
we have the following relation for the quantum potential:

$$𝓠=\frac{\alpha \mathrm{}^2}{2}\left\{4\alpha \left[\varphi _{\mu \nu }\varphi _{\mu \nu }^{(c)}\right]\left[\varphi ^{\mu \nu }\varphi ^{(c)\mu \nu }\right]\right\}$$
(49)
so the field equations are:

$$^2\varphi _{\mu \nu }+\alpha ^2\mathrm{}^2\left[\varphi _{\mu \nu }\varphi _{\mu \nu }^{(c)}\right]=0;_\mu \varphi ^{\mu \nu }=0$$
(50)
writing:

$$\varphi _{\mu \nu }=\varphi _{\mu \nu }^{(c)}+f_{\mu \nu }=\frac{C_{\mu \nu }}{r}+f_{\mu \nu }$$
(51)
where $`\varphi _{\mu \nu }^{(c)}=\frac{C_{\mu \nu }}{r}`$ with $`C_{\mu \nu }`$ constants, is the classical solution. The field equation leads to:

$$^2f_{\mu \nu }+\alpha ^2\mathrm{}^2f_{\mu \nu }=0$$
(52)
The above equation can be solved just as in the previous section. Then one must apply the gauge condition. This would leads to some relation for constants $`C_{\mu \nu }`$. The solution is:

$$\varphi _{\mu \nu }=\left(\begin{array}{cccc}\frac{2r_s}{r}+f& \text{0}& \text{0}& \text{0}\\ \text{0}& \text{0}& \text{0}& \text{0}\\ \text{0}& \text{0}& \text{0}& \text{0}\\ \text{0}& \text{0}& \text{0}& \text{0}\end{array}\right)$$
(53)
where $`f(r,\theta ,\phi )`$ is just as in the previous section, i.e.:

$$f(r,\theta ,\phi )=\underset{lm}{}\left[a_{lm}j_l(\alpha \mathrm{}r)+b_{lm}n_l(\alpha \mathrm{}r)\right]Y_{lm}(\theta ,\phi )$$
(54)
and $`r_s`$, the only non vanishing element of $`C_{\mu \nu }`$, is the Schwarzschild radius. This solution leads to:

$$h_{00}=\frac{r_s}{r}+\frac{f}{2}$$
(55)
$$h_{0i}=0$$
(56)
$$h_{ii}=\frac{r_s}{r}+\frac{f}{2}$$
(57)
So the metric $`g_{\mu \nu }=\eta _{\mu \nu }+h_{\mu \nu }`$ is given by:

$$g_{\mu \nu }=\left(\begin{array}{cccc}1\frac{r_s}{r}+\frac{f}{2}& \text{0}& \text{0}& \text{0}\\ \text{0}& 1\frac{r_s}{r}+\frac{f}{2}& \text{0}& \text{0}\\ \text{0}& \text{0}& 1\frac{r_s}{r}+\frac{f}{2}& \text{0}\\ \text{0}& \text{0}& \text{0}& 1\frac{r_s}{r}+\frac{f}{2}\end{array}\right)$$
(58)

## IV SOME OBSERVABLE RESULTS

At this end let us look for some observable results of the static spherical symmetric solution of the previous section. It is a good idea to see that the extra quantum terms how affect the classical results of the Schwarzschild metric. We present here two effrcts.

1. Light deflection

As an example consider the light bending effect. Suppose a light ray approaches a star from infinity along the $`z`$ direction and in the $`xz`$ plane with impact parameter $`b`$. The light ray path would be given by the null geodesic relation, written in the linear field approximation:

$$dk_\mu \frac{1}{2}\left(_\mu h_\beta ^\alpha \right)k_\alpha dx^\beta =0$$
(59)
where $`k_\mu `$ is the wave number. This equation simplifies to:

$$dk_x=\frac{\omega }{2}\frac{}{x}\left(\frac{2r_s}{r}+f\right)_{x=b}dz$$
(60)
where $`\omega `$ is the frequency of the light beam.

For simplicity, we assume only the term $`a_{00}`$ contribute in $`f`$ so we setIf one wants to have a metric not dependent on the spherical angles at infinity, only the coefficients of $`Y_{00}`$ are nonzero. The behaviour of $`j_0`$ and $`n_0`$ are very similar, so for simplicity we only consider the first one. Considering the second does not change the physics.:

$$f=a_{00}j_0(\alpha \mathrm{}r)=a_{00}\frac{\mathrm{sin}\alpha \mathrm{}r}{\alpha \mathrm{}r}$$
(61)
Integrating the relation (60), one gets:

$$\theta \frac{_{\mathrm{}}^+\mathrm{}𝑑k_x}{k_z}=\frac{2r_s}{b}a_{00}\left(\mathrm{cos}\alpha \mathrm{}b\frac{\mathrm{sin}\alpha \mathrm{}b}{\alpha \mathrm{}b}\right)$$
(62)
in which $`\theta `$ is the angle of bending. The first term represents the classical light bending relation, while the second term is the quantum corrections. The dependence of the bending angle on the impact parameter is plotted in the figure (1).

As it can be seen, the quantum solution is an oscillating curve around the classical curve. The amplitude of the oscillations is scaled by $`a_{00}`$, which is an arbitrary constant. So if any oscillatory behaviour be observed in the light bending effect, it can be described by the quantum effects.
2. Tidal forces

In order to investigate how the quantum effects change the curvature of the space-time, it is a good idea to derive the tidal forces. This is a simple task, noting that the tidal forces are in fact the gradient in the gravitational field in locally freely falling frame. In figures (2) and (3) the classical and quantum fieldplots of tidal forces at some point are ploted. Also to see the effect better, in figures (4) and (5) the classical and quantum shape of a sphere falling on a spherically symmetric mass distribution is shown. Note the quantum oscillations in the shape in the last figure.

Acknowledgement The first author wishes to thank Tehran University for supporting this project under the grants provided by the research council.