format: text

text:
# Spontaneous spinning of a magnet levitating over a superconductor

## Abstract

A permanent magnet levitating over a superconductor is found to spontaneously spin, overcoming resistance to air friction. We explain the physics behind this remarkable effect.

It is well known that a permanent magnet can levitate over a superconductor due to the Meissner effect. Furthermore, for type II superconductors, levitation over a flat superconductor is stable due to flux penetration and pinning of flux lineshellm . However it appears not to be widely known that under appropriate conditions a levitating magnet has a strong tendency to start spinning: independent of initial conditions the rotation will reach a terminal velocity and persist, overcoming resistance due to air friction. After coming across this phenomenon independently we found out that this remarkable effect has been seen and discussed beforeold1 ; old2 . However its physical origin has not been completely elucidated before. The reader can see examples of this remarkable effect at the website movie .

A photograph of our experimental setup is shown in Figure 1. A $`YBCO`$ disk of diameter $`22mm`$ and thickness $`3mm`$ rests on a metal base that is submerged in liquid nitrogen in a flat glass container. A permanent $`Nd_2Fe_{14}B`$ magnet disk of radius $`R=2.5mm`$ and thickness $`h=1.2mm`$ is placed vertically over the superconductor, as shown in the figures. The magnet is product number $`641895`$ from Radio Shack Corporation. The height at which the magnet levitates over the superconductor is approximately $`0.4`$ mm. It is determined by various parameters in the experimental setup such as the thickness of the superconductor, the mass of the magnet and the strength of the magnetization, as discussed in Ref.hellm . The magnetization of the permanent magnet is perpendicular to the magnet disk surface, hence rotation of the magnet around a horizontal axis going through the disk center does not change the magnetic field.

If the magnet has a high degree of homogeneity it will start rotating almost immediately, while if it is inhomogeneous it will start rotational oscillations around an equilibrium orientation which increase in amplitude, until it ’goes over the top’ and starts rotating either clockwise or counterclockwise at random depending on initial conditions. We have found this by experimentation with 8 different magnets, some with a high degree of homogeneity and some with less. The torque that makes the magnet oscillate and spin is to lowest order proportional to the magnet’s angular velocity: if the magnet is initially at rest with zero angular velocity it is in unstable equilibrium, and any small perturbation will initiate the motion .

In agreement with previous workold1 ; old2  we find that the tendency to spinning is strongest when the magnet is situated in a region of large temperature gradient. In the setup of Fig. 1 with container height $`1.2cm`$ we find that a metal base of approximately $`5mm`$ is optimal for strongest spinning tendency. For such height the lower part of the magnet is still inside the container, while the upper part is outside; this presumably gives rise to largest temperature gradient. The terminal angular velocity for that case is approximately 4 revolutions/second. For a lower metal base the tendency to spinning and the terminal velocity is reduced, and if we use no metal base in this container, or if we use a deeper container, the tendency to spinning disappears. In such cases we would expect the entire magnet to be at a more or less uniform cold temperature due to proximity to the liquid $`N_2`$. Furthermore, we found that with the metal base and the flat container the spinning can be stopped by placing another container with liquid $`N_2`$ directly above the magnet, thus cooling the environment around the upper part of the magnet. Furthermore we have found that the tendency to spinning is reduced if the magnet is placed not vertically but at an angle; in that case, the temperature difference between top and bottom of the magnet is obviously reduced.

The observed effect is surprising because there is no obvious source for the kinetic energy of rotation that the magnet acquires, or for the torque that makes the magnet spin overcoming air friction. The effect does not depend on details of the experimental setup described above. We have observed it also with other superconducting samples of different shape, even with several small pieces of superconductor supporting the levitating magnet. The effect is also seen when several of the magnets described above are put together making a cylinder of the same radius but larger height.

It was found in reference old1  that the effect does not originate in convective currents, i.e. ’wind’, originating in the evaporating $`N_2`$, by placing the magnet inside a test tube and observing that the effect persists unchanged. We have repeated the same check. Additionally it was shown in ref. old2  that the effect persists when the air pressure is reduced. Furthermore it was conjectured in refs. old1 ; old2  that the effect depends on the thermomagnetic properties of the magnet, namely on the reversible change of magnetization with temperature. However no direct experimental evidence was presented that the effect is connected to the magnetism of the magnet or the superconductivity of the superconductor. In principle one could imagine that the effect would occur also for a non-magnetic disk in a large temperature gradient if friction with the environment can be minimized.

However we have repeated the experiment with superconductors of different thicknesses, and observed a stronger tendency for spinning and a larger terminal angular velocity under the same conditions of temperature gradient for increasing thickness. This establishes that the superconductor plays a crucial role beyond simply levitating the magnet so that it can rotate without significant friction. Furthermore we have repeated the experiment with three identical magnets forming a cylinder of the same radius and three times the original height, and found that the tendency to spinning is only slightly reduced; when we then replaced the middle magnet with another magnet that had been previously heated to a temperature high enough to completely eliminate its magnetization, it was found that the tendency to spinning was now very significantly reduced. This establishes that the magnetism of the magnet plays a crucial role beyond simply levitating the magnet and that the spinning effect is larger for larger magnetization.

We have also observed that the effect depends on the thermal conductivity of the magnet. We did not have other magnets of different thermal conductivity but otherwise similar characteristics available, however, we observed that if the magnet is wrapped with thin Al foil, the effect becomes weaker or disappears altogether. The thermal conductivity of $`Al`$ is much larger than that of $`Nd_2Fe_{14}B`$. This observation suggests that the effect will only occur if the thermal conductivity of the magnet is not too large.

Because the effect appears to persists indefinitely it does not seem possible that the energy to sustain the motion is supplied by the superconductor. Hence we conclude that it originates in heat from the environment, i.e. that the rotating magnet constitutes a simple heat engine as depicted in Figure 2. At any instant when the magnet is rotating, its upper half will be at a lower temperature than its environment, because half a period earlier it was in a region of much colder temperature. Similarly, its lower half will be at temperature higher than its environment because it was in a higher temperature region half a period ago. Hence the upper half will absorb heat from the environment and the lower half will release heat to the environment. If the heat absorbed by the upper half is more than the heat released by the lower half, the difference will be converted into work.

It remains to be explained how the conversion of heat to work depicted in Figure 2 occurs. We argue that it is a necessary consequence of the existence of a temperature gradient in the region where the magnet is and of the fact that the magnetization of the $`Nd_2Fe_{14}B`$ magnet is a decreasing function of temperature in the temperature range where the magnet is, as discussed in Refs.old1 ; old2 . Qualitatively it works as follows: the superconductor exerts a force on each element of the magnet that is an increasing function of its magnetization. The magnetization in turn is an increasing function of temperature. An element of the magnet moving down has a higher temperature than its symmetric counterpart at the same distance from the superconductor that is moving up. The unequal force exerted by the superconductor on these two magnet elements leads to a net torque in the same direction as the existing angular velocity.

To make this argument concrete, consider the expression for the energy cost of putting vortices originating in a magnetic moment $`\stackrel{}{m}`$ (parallel to the superconductor) through the superconductor, given in Ref. hellm

$$E=\frac{H_{c1}L}{\pi }\frac{m}{d}C\frac{m}{d}$$
(1)
with $`H_{c1}`$ the lower critical field of the superconductor, $`L`$ the thickness of the superconductor and $`d`$ the distance from the superconductor to the magnetic moment. Let $`(r,\phi )`$ be the radial and angular coordinates of an element of magnetic moment $`\stackrel{}{dm}`$ for which $`d=d_0+rcos\phi `$, with $`d_0`$ the distance from the superconductor to the center of the magnet disk. We measure $`\phi `$ with respect to a vertical axis $`z`$. The torque on the magnet due to the force exerted by the superconductor on this element of magnetic moment is

$$d\tau =\frac{E}{\phi }=\frac{C\times dm\times rsin\phi }{(d_0+rcos\phi )^2}$$
(2)
The element of magnetic moment is given by

$$dm=M[T(r,\phi )]hrdrd\phi $$
(3)
where $`M[T]`$ is the local magnetization of the disk, which we assume to depend on the local temperature $`T(r,\phi )`$. If the magnet is rotating and has finite thermal conductivity we will have

$$T(r,\phi )T(r,\phi )$$
(4)
with $`T`$ lower for the ’rising’ side, hence there will be a net torque on the magnet. It is plausible to assume to lowest order

$$T(r,\phi )T(r,\phi )=krsin\phi $$
(5)
with $`k`$ a constant, hence the net torque is

$$\tau =Chk\frac{M}{T}_0^R𝑑r_0^\pi 𝑑\phi \frac{r^3sin^2\phi }{(d_0+rcos\phi )^2}$$
(6)
The constant $`k`$ will depend on the angular velocity $`\omega `$ of the magnet and its thermal conductivity $`\kappa `$, and in particular will go to zero both when $`\omega 0`$ and when $`\kappa \mathrm{}`$. The torque $`\stackrel{}{\tau }`$ is a vector parallel or antiparallel to the angular velocity $`\stackrel{}{\omega }`$ depending on whether $`\frac{M}{T}>0`$ or $`\frac{M}{T}<0`$ respectively. In the second case the torque opposes the motion and together with air friction will rapidly dampen any initial oscillation or rotation; in the first case instead, the torque will act in the presence of any small initial perturbation to increase the motion and the magnet will gain angular momentum and kinetic energy. In the presence of air friction the magnet will then reach a terminal angular velocity and continue its rotation indefinitely driven by this torque.

The height $`d_0`$ at which the magnet levitates is proportional to $`L^{1/2}`$hellm , with $`L`$ the thickness of the superconductor, hence from Eqs. (6) and (1) we expect the torque to increase as $`L^{1/2}`$. This is in qualitative agreement with our experimental observations discussed earlier.

The expression Eq. (6) is valid for the case where the entire magnet has magnetization with the same $`\frac{M}{T}`$. If as described in our experiment earlier, part of the magnet has no magnetization, this part does not contribute to the torque Eq. (6) but does contribute to the moment of inertia, thus reducing the tendency to spinning in agreement with our observations discussed earlier.

Next we discuss the principle on which this ’heat engine’ is based. We consider the simplified example shown in Figure 3 of a magnetized ’pendulum’ on top of a superconductor in the presence of a temperature gradient. The magnet has two possible positions, one characterized by height $`d_1`$ and environment temperature $`T_1`$ and the other by $`d_2`$ and $`T_2`$, with $`d_1>d_2`$ and $`T_1<T_2`$. Furthermore, the magnet’s magnetic moment when it is at temperature $`T_i`$ is $`m_i`$, with $`m_2>m_1`$. In the first step the magnet (at temperature $`T_1`$) falls from $`d_1`$ to $`d_2`$ and the superconductor increases its energy by

$$ϵ_1=Cm_1(\frac{1}{d_2}\frac{1}{d_1})$$
(7)
In the second step at position $`d_2`$, the magnet comes to thermal equilibrium with its environment at lower temperature $`T_2`$ by releasing heat $`q_2`$ , its magnetization increasing to $`m_2`$. The energy of the superconductor increases by

$$ϵ_2=\frac{C}{d_2}(m_2m_1)$$
(8)
In the third step, the magnet rises again to position $`d_1`$ and the superconductor lowers its energy by

$$ϵ_3=Cm_2(\frac{1}{d_2}\frac{1}{d_1})$$
(9)
and in the fourth step the magnet absorbs heat $`q_1`$ and the magnet heats to temperature $`T_1`$, and its magnetization decreases to $`m_1`$, while the superconductor decreases its energy by

$$ϵ_2=\frac{C}{d_1}(m_2m_1)$$
(10)
Note that $`ϵ_1+ϵ_2ϵ_3ϵ_4=0`$ as required by conservation of energy for the superconductor. The difference in the energy increase of the superconductor when the magnet falls and its decrease when the magnet rises is the gain in the kinetic energy of the magnet:

$$\mathrm{\Delta }K=ϵ_3ϵ_1=C(m_2m_1)(\frac{1}{d_2}\frac{1}{d_1})$$
(11)
and conservation of energy in a cycle requires that

$$q_1q_2=\mathrm{\Delta }K$$
(12)
hence the magnet absorbs more heat than it releases, and the difference is converted into work, either increasing the kinetic or potential energy of the magnet or dissipated in air friction or both.

The process described above is an idealization but describes the essence of the mechanism giving rise to the spontaneous rotation. It is a necessary condition for the above analysis to be valid that the magnet disk does not conduct heat too fast, so that the side coming down has a higher temperature than the side going up. The thermal conductivity of $`Nd_2Fe_{14}B`$ is $`\kappa =9W/mK`$, and its specific heat is $`C=0.5\times 10^7ergs/grK`$. Using these values and the dimensions of our magnet we estimate the time it takes for heat to conduct across the magnet as approximately

$$\mathrm{\Delta }t\frac{C\mu }{\pi \kappa d}2.5sec$$
(13)
where $`\mu =0.174gr`$ is the mass of the magnet. This is much larger than the observed rotation period, hence we conclude that indeed the temperature profile on the magnet is not equilibrated as required by the mechanism discussed here. Instead, the thermal conductivity of $`Al`$ is approximately $`25`$ times larger than that of $`Nd_2Fe_{14}B`$, $`\kappa =237W/mK`$. Hence when wrapped with $`Al`$ foil the magnet would have to turn faster than $`10`$ revolutions/second to prevent thermal equilibrium throughout the magnet and allow this mechanism to operate. The resistance from air friction would however be much larger for such large angular velocity, which is presumably why the effect does not occur in that case.

In summary, we have explored the physics behind the unusual observation that a magnet levitating over a superconductor is found to spontaneously spin under the right conditions. We have explained the phenomenon as constituting a simple heat engine, where heat is absorbed from a high temperature source, a smaller amount is released to a lower temperature sink, and the difference is converted into work through the interaction between the magnet and the superconductor. Compared to other simple heat engines such as Carnot’s or Stirling’s the present one appears to be simplest, as it has only one moving part, the magnet itself, and all the steps in the cycle are combined in a single process, the rotation of the magnet. Possible applications should be explored.

###### Acknowledgements.

The authors are grateful to B. Maple for kindly providing the superconducting sample and the liquid $`N_2`$ used in the experiments.


text:
# The measurements, parameters and construction of Web proxy cache

## 1 Introduction

This paper considers approaches to the problem of optimization of proxy cache construction based on theoretical and experimental study. Until recently, caching was an optional service for users who voluntary configured their browsers to redirect request through a proxy. The Internet Service Providers interpose the caching systems in the strategic places at the organization boundaries.

In experimental research of the cache systems and construction of the mathematical models should lead to the growth of the caching effectiveness with minimal financial expenditures. A better algorithm that increases hit ratios by only several percent would be equivalent to a multiple growth in cache size. In order to find the optimal cache size Kelly and Reeves  are guided on economical methods like monetary cost of memory and bandwidth.

Our analysis is based on the model presented by Breslau et al.  and extended by Wolman et al.  to incorporate the steady-state behavior and documents’ rate of change. One difficulty in the study of Web caching is that there are many cache replacement policies and many factors affecting their performance. Wolman et al. parameterized the model using population size, population request rate, document rate of change, size of object universe, and popularity distribution of objects. Several researchers  include additional factors like object size, miss penalty, temporal locality, and long-term access frequency. Our research allows the calculation of the lower limit of the cache size corresponding to aggregated bandwidth of external links when performance mount to the effective level equal to 35%.

Such an approach can be easily generalized to describe any applications based on Zipf-like distribution such as Content Distribution Networks (CDN) , peer-to-peer systems , Internet search engines, etc.

The paper investigates how the different parameters of proxy cache influence on its performance. The methods of the measurements including the treatment of experimental data and analytical formulas for calculations are discussed. The correlations between the significant parameters are investigated and the corresponding Figures and Tables are constructed. On the basis of the new experimental data the analytical model is modernized and the addition to the cache construction and to replacement algorithm is proposed.

Special consideration is given to the modification of the key formulas introduced by Wolman at al. . In order to describe the renewal effect of Web documents in the global network the alternative model is developed. The document rate of change $`\mu (i)`$ is supposed to depend on popularity index $`i`$ so the Zipf-like distribution with the new exponent $`\alpha _R`$ describes the mentioned effects of the Wolman at al. model. As a result of the special experiment the values of exponents $`\alpha `$, $`\alpha _R`$ and of the document rates of change $`\mu _p`$, $`\mu _u`$ are calculated.

## 2 Plan of the measurements

The scheme of the Web caching could be presented in the following way: there are users which send requests to the global network and receive information through a cache system as it is shown on the Fig. 1. Some documents are requested repeatedly and therefore they should be held in the cache system.

The relative frequency  of requests to Web pages follows Zipf-like distribution . This distribution states that the relative probability of a request for the $`i`$’th most popular page is

$$\theta _i=\frac{A}{i^\alpha },$$
(1)
where $`A=\theta _1`$ is the probability of the most popular item and $`\alpha `$ is a positive exponential value less then unity.

We proponed the following test: the size of proxy cache $`S_{eff}`$ was varying. These values correspond to incoming traffic for one, two, three, and six days. All types of documents both cacheable and uncacheble are taken into account for calculating the ratio $`S_{eff}/\nu _{int}`$.

The network of Samara State Aerospace University has been chosen as an experimental field. The proxy cache of SSAU is a two-processor Linux server with SQUID proxy installed. All hierarchical links were disconnected before the experiment began. The statistics of requests for the long time $`T_{st}month`$ were collected for each point.

In order to modify the model proposed by Wolman at al. we collected the requests during the time $`T_{st}t_u`$. The $`t_u`$ is the mean lifetime of those documents, who’s popularity is $`\vartheta _i=1`$ ($`\vartheta _i=\theta _ik`$).

## 3 The original results and their processing

The primary results of experiments are summarized in the Table 1, where the following notations and abridgements are used:

* The variables $`\nu _{int}`$ and $`\nu _{out}`$ are the incoming and the outgoing request’s streams of cache proxy as it is shown on the Fig. 1. They are measured as the number of users requests per day (Rpd - request per day). It should be noted that the variable $`\nu _{out}=\lambda N`$ describes the request stream from a collective user that is a significant parameter of the Wolman at al. model.
* $`H`$ is the performance of the cache system or, in other words, hit ratio.
* The same variables with upper index $`B`$ describe the system in the units of transmitting traffic ($`Kbps`$ \- Kbit per second).
* $`E(S)`$ is the mean size of documents received from the global network directly.
* $`E(C)`$ is the mean size of documents from the cache.
* Finally, the time $`T_{st}`$ corresponds to the quantity of days when statistics are collected.

The statistics collected were processed by scripts specially written for the task. Originally, so-called cacheable documents that can be stored in a proxy are selected from the general list. Later the corresponding Zipf-like distribution was constructed where the documents were placed in the order of reducing of popularity index $`\vartheta _i`$. The fragment of this list is shown bellow:

| 112 | http://www.ixbt.com/images/empty.gif | (line 457), i.e. $`\vartheta _{457}=112`$ |
| --- | --- | --- |
| 111 | http://cacheserver.myecom.net/main/images/adlogo.jpg | (line 458) |
| . | . | . |
| 2 | http://zzz.net.ru/images/spacer.gif | (line 78166=$`M`$) |
| . | . | . |
| 1 | http://www.muz-tv.ru/chat/chat-top.html | (line 200045) |

The number of unique cacheable documents or the quantity of the lines in the list mentioned above is $`p`$. The line number of the last document, which was requested twice ($`\vartheta _M=2`$), is $`M`$.

In order to find the general number of cacheable documents $`k`$ the following sum was calculated

$$k=\underset{i=1}{\overset{p}{}}\vartheta _i$$
(2)
where $`\vartheta _i`$ is the number of cache requests for the $`i`$’th most popular document. The portion of cacheable documents $`p_c`$ is defined as

$$k=p_c\nu _{out}T_{st},$$
(3)
where $`K=\nu _{out}T_{st}`$ is the general number of all documents received from the global network including both cacheable and uncacheable ones.

The value of $`\alpha `$’s was calculated using the equation

$$\alpha =12M/\underset{i=1}{\overset{M}{}}\vartheta _i=12M/(kp+M)$$
(4)
Finally, analyzing the log files mean lifetime $`t_u`$ was calculated for those documents, with a popularity of $`\vartheta _i=1`$. Such statistics also determine the lifetime $`T_{eff}`$ of cache objects with the citing index $`\vartheta _i=2`$, i.e. those items, which have been stored in proxy cache, requested one time from a proxy, and deleted subsequently (see Tab. 2).

## 4 Basic correlation

The first family of the curves, which should be analyzed, is the dependence of $`t_u`$ and $`T_{eff}`$ on cache size $`S_{eff}/\nu _{int}`$. As it was discovered in the Ref. , these curves define the ratios between elements of cache construction:

* A kernel $`S_k`$ that contains popular documents with $`\vartheta _i2`$.
* An accessory part $`S_u`$ that keeps unpopular documents requested from the Internet once, i.e. $`\vartheta _i=1`$.
* A managing part $`S_m`$ that contains statistics of requests and rules for replacement of cache objects.

Earlier  we have got the following Egs.

$`S_k`$ $`=`$ $`{\displaystyle \frac{(1\alpha )H}{2}}\nu _{out}T_{eff}`$ (5)
$`{\displaystyle \frac{S_k}{S_u}}`$ $`=`$ $`{\displaystyle \frac{T_{eff}}{(2^{1/\alpha }1)t_u}}={\displaystyle \frac{M}{pM}}{\displaystyle \frac{T_{eff}}{t_u}}`$ (6)
Analysis of the experimental data shown in the Tab. 2 and on the Fig. 2 leads to the facts that variables $`t_u`$ and $`T_{eff}`$ are directly proportional to the cache size $`S_{eff}/\nu _{int}`$ and can be considered as coincided values:

$$T_{eff}t_u$$
(7)
The only deflection indicated for effectiveness of replacement algorithms is revealed at small cache size when performance is far from the optimal value.

In other words the kernel and accessory parts are approximately correlated as 1:2 or less then 40% of storage capacity has been used for the basic goal to store the repeatedly requested documents. It is remarkable that the number of documents $`M`$, which must be stored in the cache system for one and two months traffic, is less then three and six-days incoming stream correspondingly.

The second family of curves intended for study is the dependence of the system performance $`H`$ and $`H^B`$ on its relative size $`S_{eff}/\nu _{int}`$, see Fig 3.

The hit rate $`H`$ of the web cache is considered to grow in a log-like fashion as a function of cache size . From the expression for cache performance

$$H=p_c\underset{1}{\overset{S_k}{}}\frac{A}{x^\alpha }𝑑x,$$
(8)
the following dependence appears

$$\frac{H_1}{H_2}=\left(\frac{S_1}{S_2}\right)^{1\alpha },$$
(9)
that allows us to talk about power fashion.

The Fig 3 illustrates the fact that the dependence between $`H`$ and $`S_{eff}/\nu _{int}`$ is successfully described by Eq. (9) with $`\alpha =0.77`$ and the curve $`H^B`$ is not predictable. This effect needs additional study especially because the Tab. 1 shows a positive difference between the mean size of cacheable documents $`E(S)`$ and the mean size of all items $`E(C)`$.

## 5 Renewal of Web documents

Rapid development of computer technology at the end of the last millennium lead to the appearance of a virtual world with its own laws. Unfortunately, during this period little attention was given to studying of fundamental principles of virtual life.

Nowadays we can afford some respite and some resources should be transferred to studying of the fundamental laws and to optimization of the vital systems on the basis of recent knowledge. Frankly speaking, any research area may be considered as a scientific field if and only if its basic principles are enveloped in mathematical form and new rules can be predicted on the basis of confirmed facts.

Our analysis is based on the model presented by Breslay et al  and extended by Wolman et al.  to incorporate document rate of change. The Wolman model yields formulas to predict steady-state properties of Web caching systems parameterized by population size $`N`$, population request rates $`\lambda `$, document rate of change $`\mu `$, size of object universe $`n`$ and a popularity distribution for objects.

The key principles of their theoretical study has been developed in our work  to describe the basic effects as a ratio between construction parts, steady-state performance, optimal size, etc. in an alternative way. We want to outline that the results received and the plan of the experiment follow from the theoretical model .

The key formulas from Wolman et al.

$`C_N`$ $`=`$ $`{\displaystyle \underset{1}{\overset{n}{}}}{\displaystyle \frac{1}{Cx^\alpha }}\left({\displaystyle \frac{1}{1+\frac{\mu Cx^\alpha }{\lambda N}}}\right)𝑑x`$ (10)
$`C`$ $`=`$ $`{\displaystyle \underset{1}{\overset{n}{}}}{\displaystyle \frac{1}{x^\alpha }}𝑑x`$ (11)
yield $`C_N`$, the aggregate object hit ratio considering only cacheable objects. Document rates of change $`\mu `$ were considered to take two different values, one for popular documents $`\mu _p`$ and another for unpopular documents $`\mu _u`$. An additional multiplier from Eq. (10) throws off one effective query to any document from cache during time $`T_{ch}=1/\mu `$ between its changes. Such an updating request must be redirected to the global network to get the renewed Web page.

In this paper we assume that the document rate of change $`\mu (i)`$ depends on its popularity $`i`$. The mathematical equivalent of this assertion is that the steady-state process is again described by Zipf-like distribution with a small $`\alpha _R`$ as it is shown of the Fig. 4.

The difference between an ideal performance of a cache system and its real value has been conditioned by the renewal of documents in the global network

$$\mathrm{}H=\left(\underset{i=1}{\overset{M}{}}\vartheta _iHK\right)/K$$
(12)
and can be investigated in an experimental way. Here

$$\mathrm{}k=\underset{i=1}{\overset{M}{}}\vartheta _iHK=kHK+Mp,$$
(13)
$`\mathrm{}k`$ is the number of updating requests conditioned by renewal of Web documents and

$$k_R=HK+pM$$
(14)
Therefore we have to explore the log file collected for the time $`T_{st}T_{eff}\mathrm{}T`$ using the data from the Tab. 2. Then the condition $`M(T_{st})S_k`$ has been fulfilled and all repeated requests are explained by the renewal effect. The ideal Zipf-like distribution corresponds to the upper line on the Fig. 4 with exponent $`\alpha `$ from Eq. (4), but the real hit ratio $`H`$ determines $`\alpha _R`$ as

$$\alpha _R=12M/HK$$
(15)
Now it is easy to find $`\mu _i`$

$$\mu (i)=\frac{\vartheta _i^{id}\vartheta _i^R}{T_{st}}=\frac{(p/i)^\alpha (p/i)^{\alpha _R}}{T_{st}},$$
(16)
or

$$\mu (i)=\frac{1}{T_{st}}\frac{1(i/p)^\mathrm{}\alpha }{(i/p)^\alpha }.$$
(17)
We can consider $`\mu _u=\mu (p/4)`$ and $`\mu _p=\mu (p/100)`$, then the results can be summarized in the Table 3.

The key difference between our analytical model and presented by Wolman at al. model are:

* The document rate of change $`\mu (i)`$ is described by continuous variables that depend on popularity $`i`$. Wolman et al. assume that the document rate of change $`\mu `$ takes two different values, one for popular documents $`\mu _p`$ and another for unpopular documents $`\mu _u`$.
* The renewal of Web documents is incorporated in a steady-state process as a Zipf-like distribution with a small $`\alpha _R`$.
* We assume also that size of cache system is restricted.
* A considerable part of cacheable documents $`pM`$ are requested from the global network only one time accordingly to Zipf-like distribution. The more accurate expression for an ideal hit ratio $`H_i`$ has been shown in Ref.

$$H_i2^{(\alpha 1)/\alpha },$$
(18)
It has to be modified to taking into account the renewal effect:

$$H_i2^{(\alpha 1)/\alpha }(1\alpha )/(1\alpha _R),$$
(19)

Two types of correlations pretend to a role of fundamental laws that describe cache systems :

* A Zipf-like distribution, see Eq. (1)
* Normalizing conditions or a sum of the probability to request the universe of $`1nk`$ objects.

The above mentioned laws could be applied to the special points of Zipf distribution and two of them, $`M`$ and $`p`$, are used for the construction of the theory. Then the Zipf-like distribution leads to

$`{\displaystyle \frac{Ak_R}{M^{\alpha _R}}}`$ $`=`$ $`2,`$ (20)
$`{\displaystyle \frac{Ak_R}{p^{\alpha _R}}}`$ $`=`$ $`1.`$ (21)
Normalizing conditions for the first $`M`$ and $`p`$ documents from cache give

$`{\displaystyle \underset{1}{\overset{M}{}}}{\displaystyle \frac{A}{x^{\alpha _R}}}𝑑x`$ $`=`$ $`H_i`$ (22)
$`{\displaystyle \underset{1}{\overset{p}{}}}{\displaystyle \frac{A}{x^{\alpha _R}}}𝑑x`$ $`=`$ $`1`$ (23)
Here $`H_i`$ is an ideal (steady-state) performance for cacheable documents. For a real system the Eq. (22) has been transformed to

$$H=p_c\underset{1}{\overset{S_k}{}}\frac{A}{x^{\alpha _R}}𝑑x,$$
(24)
where $`S_k`$ is the number of cache objects in the kernel.

It is important that Zipf exponent $`\alpha `$ grows as the experiment time $`T_{st}`$ increases.

$$\mathrm{}\alpha =\alpha (T_{st}^2)\alpha (T_{st}^1)=C\mathrm{ln}(T_{st}^2/T_{st}^1)$$
(25)
## 6 Summary and future work

The aim of this paper is a study of a Web cache system in order to optimize proxy cache systems and to modernize construction principles. Our investigations lead to the criteria for the optimal usage of storage capacity and to allow to description of the basic effects as the ratio between the construction parts, steady-state performance, optimal size, etc. We want to outline that the results received and the plan of the experiment follow from the theoretical model.

Special consideration is given to the modification of the key formulas supposed by Wolman at al. . The document rate of change $`\mu (i)`$ is supposed to depend on popularity index $`i`$ so the Zipf-like distribution with the new exponent $`\alpha _R`$ describes the effects of the renewal of Web documents.

The main result of any research of cache system is finding a way for increasing hit ratio. We can conclude that with growth of $`S_k/S_u`$ the hit rate $`H`$ increases. A general feature of the current algorithm is the only documents requested two times and more during the time $`t_u`$ are included in the kernel $`S_k`$. It is a fact that $`M`$ less then the cache size $`S_{eff}`$, i.e. all necessary items could be stored in a cache system at the second, third and forth experimental points on the Tab. 2. Therefore we can make the conclusion that a caching algorithm based on a rigid tie of the mean parameters to the time $`t_u`$ is ineffective.

One method of resolving the current situation is reconstructing of cache systems. Such a construction suggested in the Sec. 4 must be implemented to provide a rigid ratio between its elements. The requests’ statistics must be kept even for those cacheable documents that have been deleted from the cache. This statistics should exist for long time. This is one of the key differences between our construction and the existing ones that usually operate only with the current set of documents in cache. An inseparable part of the new construction is the replacement algorithm based on Zipf law. Now we have prepared the corresponding application for an international patent.

An additional question, which will demanded a new experimental search, is the dependence that describes the renewal of Web documents. Probably, the ratio $`\mathrm{}\alpha /\alpha `$ could be considered as a universal constant. The next direction for our development plans is an investigation of a model of cache interaction in hierarchical caching system.


text:
# Molecular conduction: paradigms and possibilities

## I Molecular conduction: what is the underlying physics?

Recently several researchers have measured charge transport in single or small groups of organic molecules connected to metal contacts . In parallel, there have been theoretical attempts at understanding molecular conduction, both at the semi-empirical  and first-principles  levels. Understanding molecular conduction is challenging, since it involves not just the intrinsic chemistry of the molecule, but extrinsic factors as well, such as the metal-molecule bonding geometry, contact surface microstructure and the electrostatics of the environment. The aim of this article is to discuss the various physical factors that influence molecular current-voltage (I-V) characteristics, and our attempts to model them both qualitatively and quantitatively.

A typical two-terminal molecular I-V looks like Fig. 1, often with a clear conductance gap . How does one understand such an I-V? The first step is to draw an energy-level diagram, as in Fig. 2. An isolated molecule has a discrete set of energy levels, with a highest occupied (HOMO) and a lowest unoccupied molecular orbital (LUMO), separated by a HOMO-LUMO gap (HLG). On connecting the molecule to metallic contacts, two changes happen: (i) the discrete molecular levels broaden into a quasicontinuum density of states (DOS) due to hybridization with the metal wave functions. Often the DOS retains a distinct peak structure, in which case it is still useful to think in terms of broadened molecular energy “levels”; (ii) the difference in work functions between the molecule and the metal causes charge transfer and band alignment between the two materials. The molecule equilibrates with the contact with an overall chemical potential set by the metal Fermi energy $`E_F`$, typically lying inside the HLG. Under an applied bias the molecule tries to equilibrate simultaneously with both contacts with bias-separated chemical potentials $`\mu _{1,2}`$, and is thereby driven strongly out of equilibrium. As long as the bias is small and $`\mu _{1,2}`$ lie in the HLG, the HOMO levels stay filled and LUMOs are empty and there is no current. However, when the bias is large enough that either $`\mu _1`$ or $`\mu _2`$ crosses a molecular level $`E_{\mathrm{MOL}}`$, that level is filled (reduced) by one contact and emptied (oxidized) by the other (Fig. 2a) and therefore starts conducting current . For opposite bias, the same level starts conducting when crossed by the other contact chemical potential (Fig. 2b). The net result is that for a spatially symmetric molecule with symmetrically coupled contacts the total conductance gap is given by $`4\left(E_FE_{\mathrm{MOL}}\right)`$. Molecular conduction thus depends on both the intrinsic molecular chemistry through $`E_{\mathrm{MOL}}`$ and the contact microstructure through $`E_F`$.

The intrinsic chemistry of an isolated molecule can be handled with sophisticated quantum chemical codes that can be purchased or even downloaded from the Internet. Given an appropriate basis-set (for instance, a minimal STO-3G basis) and an appropriate model for electron-electron interactions (based on first-principles density functional, Hartree-Fock or semi-empirical Hückel methods), such a code starts with an initial guess density matrix $`\rho `$ to obtain a Fock matrix $`F`$ (Fig. 3a). It then fills up the corresponding energy eigenstates with a given number $`N`$ of electrons according to equilibrium statistical mechanics to reevaluate $`\rho `$, recalculates $`F`$ and so on, until self-consistent convergence. Our molecular system differs in two ways: (i) it is open, with a varying, fractional occupancy of electronic levels; (ii) it is trying to equilibrate under bias with two different contact chemical potentials and is therefore driven strongly out of equilibrium. To solve this problem, we modify the above self-consistent scheme, as shown in Fig. 3b. The initial step, solving for the Fock matrix $`F`$, is kept unchanged from the usual prescriptions in molecular chemistry. In this step, one can use semi-empirical tight-binding/Hückel-based methods , or exploit the sophisticated numerical prowess of a standard quantum chemical package such as GAUSSIAN’98  to employ a density functional theory (DFT)-based method for evaluating an $`F`$ matrix. The $`F`$ matrix is then supplemented with self-energy matrices $`\mathrm{\Sigma }_{1,2}`$ describing an open system connected to the two contacts and involving the detailed contact microstructure, while the nonequilibrium (transport) part is set up using the nonequilibrium Green’s function (NEGF) formalism . We refer the reader to our past work for details of the NEGF equations and the calculation of the self-energy matrices . In this article, we will concentrate mainly on the physical insights.

The coupled DFT-NEGF formulation of molecular electronics is in effect the generalization of the coupled Poisson-hydrodynamic equations used extensively in analyzing device transport . We supplement the Poisson (Hartree) term with exchange-correlation corrections that are small in macroscopic devices but significant in molecules, while the semiclassical hydrodynamic equation is replaced with a fully quantum mechanical NEGF formalism .

The self-consistent formalism described above is completely general, and can be employed in modeling transport through a wide variety of physically different systems. For instance, we have used this scheme to obtain semi-empirical  and first-principles density functional (DFT)-based (LANL2DZ/B3PW91)  I-V characteristics of metallic quantum point contacts as well as semi-conducting aromatic thiol molecules bonded to Au(111) surfaces. Whether one is dealing with a carbon nanotube, a ballistic MOSFET , a spin transistor, a resonant tunneling diode  or a molecular wire, the above formalism holds. One needs simply to evaluate each of the following quantities: (i) an appropriate Fock matrix $`F`$ describing the device; (ii) a contact Fermi energy $`E_F`$ relative to which the device levels are evaluated and which determines the electrochemical potentials $`\mu _{1,2}`$; (iii) a self-consistent potential $`U_{\mathrm{SCF}}`$ describing charging effects and the electrostatic influence of the environment (this term is incorporated into the effective $`F`$ matrix of the device), and (iv) a set of self-energy matrices $`\mathrm{\Sigma }_{1,2}`$ that describe the coupling of the device with the contact, the matrices depending on the contact surface Green’s function as well as the device-contact bonding geometry. Additional self-energy matrices can be introduced to describe scattering, by phonons or polarons for example. Within the same self-consistency scheme and NEGF prescription, one can then get qualitatively different I-V characteristics just by varying $`F`$, $`E_F`$, $`U_{\mathrm{SCF}}`$ and $`\mathrm{\Sigma }`$.

The effect of each of the above parameters on the molecular I-V is schematically shown in Fig. 1. The conductance gap depends on $`\left(E_FE_{\mathrm{MOL}}\right)`$, and the maximum current level is set by the parallel combination $`\mathrm{\Gamma }_{\mathrm{eff}}=\mathrm{\Gamma }_1\mathrm{\Gamma }_2/(\mathrm{\Gamma }_1+\mathrm{\Gamma }_2)`$ of the individual contact broadenings $`\mathrm{\Gamma }_{1,2}=i(\mathrm{\Sigma }_{1,2}\mathrm{\Sigma }_{1,2}^{})`$. The current rises over a width which depends on the total broadening of the molecular levels, which in turn depends on (a) the thermal broadening $`k_BT`$, (b) the series combination $`\mathrm{\Gamma }=\mathrm{\Gamma }_1+\mathrm{\Gamma }_2`$ of the contact broadenings $`\mathrm{\Gamma }_{1,2}`$ and (c) the Coulomb charging energy $`U_0`$ to add an extra charge to the molecule (the charging energy tends to drag out the conductance peak, so that for appreciable charging energies, the contact chemical potentials may not be able to cross a molecular level easily under bias, resulting in a relatively featureless I-V characteristic). Finally, the potential profile across the molecule sets the overall voltage division factor $`\eta `$  which determines the prefactor in the ratio between the conductance gap and $`E_FE_{\mathrm{MOL}}`$ ($`\eta =0.5`$ and the prefactor equals 4 for symmetric coupling, as in Fig. 1). The voltage division factor $`\eta `$ depends on the contact geometries and characterizes the Laplace part of the potential, while the Poisson part describes self-consistent charging effects, and is characterized by $`U_0`$.

We summarize below the most challenging and physically relevant questions in obtaining a molecular I-V characteristic:

* Where is the contact Fermi energy relative to the molecular levels? $`(E_F,E_{\mathrm{MOL}})`$
* What is the broadening due to the contacts? $`(\mathrm{\Gamma }_{1,2})`$
* What is the spatial profile of the Laplace potential? $`(\eta )`$
* What is the charging energy? $`(U_0)`$

Note that each of the above quantities is in general a complicated matrix that can be modeled independently using either semi-empirical or first-principles methods. However, in order to develop a “feel” for how these quantities affect the molecular I-V as in Fig. 1, we will try to capture their essence in terms of a few characteristic scalar parameters, as defined above. We will now address the influence of each parameter one at a time below.

## II Molecular conduction: how can we model it?

### A Where is the contact Fermi energy relative to the molecular levels?

This is probably the most challenging problem to sort out. One needs to start by modeling the quantum chemistry of the molecule. For our candidate molecule phenyl dithiol (PDT), shown in Fig. 4, it is believed that the sulphur atoms bond with a Au(111) surface by desorption of the end hydrogen atoms that are then replaced by a triangle of gold atoms with sulphur sitting above their centroid . The energy levels of the isolated PDT molecule compare well with those obtained by replacing each H-atom by three gold atoms (the gold atoms introduce some additional localized levels in the HLG). Replacing the gold cluster with a self-energy describing metallic gold broadens the molecular levels into a quasicontinuous spectrum, with the localized levels developing into metal-induced gap states (MIGS) decaying spatially away from the contacts into the molecular center. Now, the energy levels of the isolated molecule itself depend sensitively on the method of calculation (a comparison plot is shown in Fig. 4). Different theoretical groups have adopted different first-principles schemes in their analysis , so the unanswered question at this point is: which method is appropriate for calculating the single-particle energy levels of an open subsystem under bias?

What is consistent among the various methods of calculation, including semi-empirical (Hückel-based) theories, is the overall shape of the orbital wave functions; for instance, the HOMO is largely sulphur-based and delocalized while the LUMO is ring-based and localized. This would seem to suggest a broad HOMO DOS and a sharp LUMO DOS, although theories don’t seem to agree even on this point (see for example Fig. 22 in  and Fig. 3 in ). The difference could arise due to different bonding geometries assumed at the gold-thiol interface. Thus the intrinsic chemistry of the molecular system needs to be cleared up, and theoretical agreement reached on the equilibrium molecular bonding properties before the transport-related issues can be sorted out.

The position of the contact Fermi energy $`E_F`$ relative to the molecular energy levels is also an unsettled issue. The Fermi energy depends sensitively on the specific model for the contact geometry. Different models for the contact (for instance DFT-based , Bethe lattice , jellium ) can give different charge contents and level broadenings of the molecule, as well as different work functions for the bulk electrodes. This could cause an appreciable shift in $`E_F`$, given the rather small gap density of states in the molecular HLG. Given that for extensively studied systems such as metal-semiconductor interfaces the precise location of $`E_F`$ is still an active topic of research , perhaps the best one can do at this point is to inquire if $`E_F`$ is closer to the molecular HOMO or LUMO level, the analogous question for a semiconductor being whether it is p-type or n-type.

Conceptually the cleanest way to address the equilibrium Fermi energy problem is by including a few layers of the contact as a cluster in a “supermolecule”, which would act as the device under investigation. Thereafter the contact is assumed unaltered during conduction, with all the “action” lying in the device sector. The advantages of including such a cluster are enormous (for a discussion, see ), such as the automatic inclusion of image charges (the supermolecule is charge neutral), avoiding uniqueness issues related to partitioning in a non-orthogonal, non-localized, atomic basis set and the proper treatment of the surface physics. Ideally, the cluster size should be significantly larger than the atomic Debye length of the contact material, while for practical purposes, it is usually limited by computational resources. $`E_F`$ is usually set by the HOMO of the (large) supermolecule, while the molecular levels can be identified by either plotting the wavefunctions or by computing the local density of states (LDOS) on the molecule. To employ this scheme to sort out the molecular chemistry and energy level structure, it is essential that the contact cluster and the molecule be calculated using the same scheme (DFT/tight-binding, etc). Attempts at performing such a computation at the semi-empirical  and DFT  levels have yielded a Fermi energy quite close to the HOMO level for PDT-Au(111) heterostructures.

The experimental situation is somewhat unclear. The conductance gap for PDT itself is different for different experimental geometries. The gap is around 3 volts for breakjunction measurements by Reed et al.  and around 4 volts for STM measurements by Hong et al. , while corresponding breakjunction measurements by David Janes’ group at Purdue indicate featureless I-V characteristics with no discernible conductance gap . Since the gap depends on $`E_F`$ which is likely to be different for the two experimental contact geometries, and could further be compromised by the presence of charging in the system, such a difference is not surprising. It is also not clear whether the conduction is through a single molecule bridging the junctions, or a series combination of molecules attached separately to the two junctions . It seems sensible therefore to treat $`E_F`$ as a fitting parameter, in the absence of precise characterization of the contact surfaces and molecular geometry. Alternatively, one could dictate the position of $`E_F`$ relative to the levels, guided by separate equilibrium cluster calculations, as discussed earlier.

Experiments incorporating a third terminal (gate) can help clarify some of these issues appreciably. For instance, a positive gate voltage lowers the molecular levels so that the Fermi energy approaches the LUMO and moves away from the HOMO. For a purely electrostatic gate control mechanism, a measured decrease in conductance will suggest HOMO (p-type) conduction, while the reverse result indicates LUMO (n-type) conduction.

Summary: Need to model the molecule and the contact bonding self-consistently within the same scheme, doing justice to the molecular quantum chemistry as well as the contact surface microstructure. The method of calculating the energy levels or the Fermi energy is still an unresolved issue.

### B What is the broadening due to the contacts?

Although experimental knowledge of the contact conditions is difficult to access, could one at least hope to model a particular idealized contact geometry and obtain an appropriate self-energy? We obtain the self-energy matrices $`\mathrm{\Sigma }_{1,2}`$ formally by an exact partitioning of the infinite metal-molecule-metal system, projecting its single particle Green’s function onto the device subspace . $`\mathrm{\Sigma }_{1,2}`$ depend on the contact surface Green’s function and the contact-molecular bondings. We obtain the couplings at the surface by simulating a large cluster from the contact coupled to the device and calculating its overlap and Fock matrices. The contact surface Green’s function is calculated using a recursive Green’s function technique taking the full group theory of the FCC Au(111) crystal surface into account .

One can replace the partitioning scheme with a scattering-formalism  that deals with the entire infinite system. Ideally, both methods (scattering formalism and the NEGF prescription) should yield the same answer; however, there is a conceptual simplification in partitioning the problem into a “device” part involving the electronically active molecule, and a “contact” part determining the lead-molecular interactions. These involve two entirely different areas of research, quantum chemistry and surface physics, so partitioning allows us to improve modeling each of them independently. Moreover, NEGF naturally allows us to include incoherent processes, which can be important even for a short molecule if there are localized states that cannot be populated from the contacts .

How do we know if we have modeled our contacts correctly? The overall “shape” of the molecular I-V can be obtained approximately without getting the bonding or the quantum chemistry right. One excellent benchmark is the quantum point contact (QPC) , the I-V characteristic of which is experimentally measured to be ohmic, with a conductance quantized in units of $`G_0=2e^2/\mathrm{}77\mu `$S . Starting with a six atom gold chain coupled to Au(111) contacts describing a QPC, we get a conductance quantized ohmic I-V, as expected . This is highly nontrivial, because conductance quantization arises out of a molecule that is perfectly transmitting over a band of energies between $`\mu _1`$ and $`\mu _2`$ (aside from Fabry-Perot type oscillations). This requires the self-energy matrices to couple the wire with the contacts seamlessly without introducing any spurious reflections. To illustrate the sensitivity of the quantization on the coupling, we scaled the overall matrix elements by a factor of five; the resultant I-V ceases to be ohmic, and resembles that of a resonantly conducting system such as PDT (Fig. 5 (i)).

The above exercise is a good check of the accuracy of the contact broadenings. Once the surface Green’s functions are deemed to be correct, it is an easy matter to replace the Au6 molecular cluster with the molecule of choice and proceed with calculating its I-V.

Summary:The QPC can be used as a benchmark for testing out the self-energy matrices. The couplings at the surface and the contact surface Green’s functions need to be calculated accurately, including the overall group theory of the metal crystal comprising the leads.

### C What is the spatial profile of the Laplace potential?

The electrostatic potential profile across the molecule can be separated into two parts: (i) the Laplace part describes the influence of the contact geometries in the absence of charges and screening effects; (ii) the Poisson part involves screening by the charges, and is determined by the charging energy of the molecule. In this section, we will concentrate on the Laplace part, and address the charging-related issues associated with the Poisson part in the next section.

The Laplace part of the potential profile is set by the relative capacitances of the contacts, and can be determined rigorously by solving the 3-D Laplace equation with the correct potential boundary conditions for the contacts. The Laplace part is characterized by the voltage-division factor $`\eta `$  which describes the proportion in which the applied voltage drops across the various contact-molecular interfaces. A convenient way to analyze the potential profile across the molecule is to incorporate the voltage-division factors in the definitions of the contact electrochemical potentials $`\mu _{1,2}`$, which would keep the molecular levels themselves fixed under drain bias in the absence of charging (next section) and move them under gate bias alone. For a two-terminal device appreciable current flow requires the capacitive couplings with the source and drain electrodes to be roughly equal (although their resistive (quantum) couplings $`\mathrm{\Gamma }_{1,2}`$ could still be quite different), leading usually to $`\eta 0.5`$.

Although the Laplace part involves essentially nineteenth century classical electrostatics, it can substantially influence device I-V characteristics. In a gated three-terminal device, for instance, a good gate control mechanism in a well-designed ballistic conductor essentially involves trying to keep the charge density near the source end of the conductor constant by pinning the device DOS to the source chemical potential . This gives an effective $`\eta 1`$, so that as the drain chemical potential ventures into the HLG under bias, the current starts to saturate . Note that this saturation mechanism is entirely different from saturation in two-terminal molecular devices, which occurs when either contact chemical potential has just crossed a molecular level and another level has not yet kicked in. In contrast to the two-terminal I-V, the three-terminal characteristic is asymmetric with respect to source-drain bias. Since the gate determines the position of the equilibrium Fermi energy through the Laplace solution, it leads to gate-modulation of the ON current in such devices.

A gate can influence the molecular electronic properties in a variety of ways: it could have a purely electrostatic effect on the channel charge and correspondingly the molecular levels, as described above. However, as pointed out by Damle et al. , good electrostatic gate control is not possible for a 10 $`\AA `$ molecule such as PDT, unless the gate oxide is prohibitively thin. Additionally, the gate can alter the properties of the contact-molecular interfaces (Schottky-barrier type effects ), or even alter the conformations of the active molecule , all of which could affect the shape of the molecular I-V.

Another example where the Laplace solution itself can influence the molecular I-V characteristic involves conduction mechanisms that require the alignment/misalignment of energy levels localized on different parts of the molecule. For instance, the Aviram-Ratner mechanism  involves a donor-bridge-acceptor system where the Laplace part of the self-consistent potential aligns the levels at the two ends for positive bias, and misaligns them for negative bias, leading to a strongly rectifying I-V characteristic. A similar example involves a quantum point contact (QPC) with a stretched bond in the middle (Fig. 5(iii)). The defect disconnects the LDOS on its two sides, allowing them to separately equilibrate with the two contacts and follow their respective chemical potentials under source-drain bias. Within the window set by $`\mu _1`$ and $`\mu _2`$, the Laplace solution causes the LDOS on the two sides to slide past each other (figure inset). Since this amounts to two transmission peaks sliding in and out of resonance, the resulting I-V shows a weak negative differential resistance .

Summary: The Laplace part of the 3-D electrostatic potential profile needs to be calculated using the boundary conditions set by the electrodes. The Laplace potential can significantly affect the molecular I-V characteristic, by aligning or misaligning different parts of the molecular LDOS or by pinning the molecular DOS to one of the electrodes in the presence of a gate terminal.

### D What is the charging energy?

The Poisson part of the potential profile carries information about screening and charging inside the molecule, and is characterized by $`U_0D_0=C_Q/C_\mathrm{\Sigma }`$, where $`U_0=e^2/C_\mathrm{\Sigma }`$ is the Coulomb charging energy, $`D_0`$ is the molecular gap DOS, the quantum capacitance $`C_Q=e^2D_0`$, and $`C_\mathrm{\Sigma }`$ describes the total electrostatic capacitive coupling to the various electrodes. The charging energy $`U_0`$ describes the “ease” with which the molecular levels can be filled or emptied, and tends to drag out the I-V characteristic as shown in Fig. 1. Note that the effect of charging can be looked upon as a voltage-dependent $`\eta (V)`$. The net capacitance $`C_\mathrm{\Sigma }`$ contributing to $`U_0`$ is determined by the geometry and dielectric constant of the molecule and the electrodes. For instance, while 10 $`\AA `$ InAs spherical quantum dots have low charging energies $`100`$ meV , a 10 $`\AA `$ isolated PDT molecular wire has a much larger charging energy $`34`$ eV, making it a lot harder to cross a level with a drain bias. At even higher charging energies $`U_0\mathrm{\Gamma }_{1,2}`$, one can get many-body effects such as Coulomb Blockade and Kondo resonance. While some Coulomb-Blockade type effects can be captured within an effective one-particle self-consistent field model within an unrestricted calculation  (such as unrestricted Hartree-Fock or spin density functional theory), doing justice to these problems requires us to go beyond the one-particle prescription we have used so far.

The Poisson solution describes the effect of adding or removing charge from the molecule, as well as the effect of redistribution of charge within the molecule, responsible for screening of the applied voltage. The efficiency of the screening process depends on the amount of material available for the reorganization of the charges. For a molecular wire much thinner than the Debye length, the net electrostatic potential is essentially given by the Laplace solution of the previous section, leading to a ramp-like potential profile . In contrast, a thick metallic wire allows sufficient screening, yielding a potential profile that is essentially flat . Such a flat potential profile can be obtained even with a thin molecular wire if the latter is embedded in a dielectric medium as part of a self-assembled monolayer (SAM), in which case the neighboring wires screen the potential profile. One way to model a SAM would be disallow any transverse variations in charge or electrostatic potential, which would amount in effect to solving the 1-D Poisson equation. This would give a highly screened potential profile even with a thin molecular wire .

Charging can lead to very interesting effects, such as the creation of an asymmetric I-V with a spatially symmetric molecule . Consider a symmetric molecule with unequal resistive (quantum) couplings to the two contacts ($`\mathrm{\Gamma }_1\mathrm{\Gamma }_2`$). Near the onset of current conduction through a HOMO level, a negative bias on the strong contact keeps it filled, while a positive bias empties it. Since the molecule gets positively charged one way but not the other, the I-V is dragged out asymmetrically such that one gets a lower current for positive bias on the stronger contact. Interestingly for conduction through a LUMO level, the sense of the I-V asymmetry reverses , because one now needs to fill the LUMO to charge up the molecule. This allows us to identify the nature of the conducting molecular orbital, which is important given that different orbitals conduct quite differently. For PDT, STM data  (the STM tip being the weaker contact) seems to indicate HOMO-based conduction. Our results (Fig. 5 (ii)) qualitatively match the I-V characteristics obtained by Reichert et al. , with an initially symmetric I-V that turns weakly asymmetric on drawing either contact away from the molecule.

The total electrostatic potential profile is the sum of the Laplace and Poisson parts, which have different effects on the molecular I-V. For typical break-junction/STM measurements on molecules such as PDT, $`C_Q/C_\mathrm{\Sigma }=U_0D_01`$, and the overall potential profile is essentially given by the Laplace potential with the net voltage dividing up as the capacitance ratio $`\eta `$ between the electrodes. If however the molecule is so strongly coupled to the substrate that there are appreciable MIGS, $`D_0`$ becomes large enough that $`C_Q/C_\mathrm{\Sigma }`$ starts to become important. In the limit of $`C_Q/C_\mathrm{\Sigma }1`$, the I-V gets dragged out substantially by charging (Fig. 1), and the final self-consistent voltage divides between the source and drain contacts according to their resistance ratio $`\mathrm{\Gamma }_2/\left(\mathrm{\Gamma }_1+\mathrm{\Gamma }_2\right)`$. This situation arises for metal nanoclusters probed by an STM tip. Since the resistance is much larger at the STM end, much of the applied voltage drops across the STM-cluster gap, so that the cluster levels remain pinned to the substrate, and the STM chemical potential alternately scans the HOMO and LUMO levels under opposite drain bias.

Summary: The charging energy can turn an otherwise symmetric I-V into an asymmetric one. Given a spatially symmetric molecule, we predict a larger current for positive bias on the stronger contact and HOMO conduction, while for LUMO conduction, the sense of the I-V asymmetry is reversed.

## III Computational issues

A brief discussion of computational issues is perhaps in order. The calculation of molecular conductance requires two steps, (a) calculating a Fock matrix given a density matrix, and (b) calculating a density matrix from the Fock matrix. The first step is the most computationally challenging part, involving the evaluation of DFT-based exchange-correlation integrals which are quite numerically complicated, especially in sophisticated basis sets involving relativistic core pseudopotentials. We find our own LDA calculation of the molecular Fock matrix in a minimal basis set to be comparable in accuracy, but not in speed, with GAUSSIAN98. Therefore we let GAUSSIAN98 do this part, exploiting decades of development that have gone into it. The second step requires contact self-energy matrices which can be calculated once and for all for realistic contact surface structures (e.g., Au(111) contacts) using a real space recursive Green’s function technique. The contact coupling matrices can be simulated in GAUSSIAN 98 with a finite-sized cluster, with care exercised to eliminate edge effects on the structure of the contact surface Green’s function  (a suitable localized basis describing gold would sort out this problem automatically).

The computational challenge for us is to solve the NEGF equations, requiring us to find the number of electrons on the molecule. Such a requirement amounts to integrating the nonequilibrium electron DOS all the way from the bottom of the contact band to the Fermi energy. Since the Green’s functions entering the DOS expression are highly peaked around the molecular levels, this process involves integrating over an energy range of several hundred volts with a millivolt accuracy for each bias point and each step of the self-consistent procedure. We have addressed this problem in two ways: (i) assuming a weak (in practice, constant) energy-dependence of the self-energy matrices (valid for Au(111) which has an energy-independent DOS near $`E_F`$), one can then perform the integrals analytically ; (ii) the nonequilibrium DOS can be divided into two groups, one lying outside the domain of $`\mu _{1,2}`$, and the rest inside that window. The first part can be integrated using a contour integration scheme , while the part between $`\mu _{1,2}`$ can be calculated either by brute force grid-based integration over the finite range between $`\mu _1`$ and $`\mu _2`$, or by reverting to the constant $`\mathrm{\Sigma }`$ approximation. Although our calculations are performed with LANL2DZ basis sets which are somewhat delocalized, it is preferable to employ relatively localized basis sets in order to avoid issues related to partitioning and to be consistent with the tight-binding approach that we are using here . Different approximations give different values of the total electron count, thereby affecting the equilibrium Fermi energy position, which clearly needs more attention. However, the approximations allow us to obtain first-principles DFT-based I-Vs for a molecule like PDT in a few hours on a SUN workstation, taking into account the details of the contact geometry.

## IV Issues we haven’t covered

We have seen that by appropriately modeling each experimental geometry, we can get qualitatively and quantitatively different conductance properties, ranging from ohmic to rectifying, switching and saturating I-V characteristics. We now outline three of the issues that we have ignored so far, namely, Conformation, Incoherence and Correlation.

(i) Conformation. One of the principal advantages of a molecule is that it is semi-flexible. This means that the conformation of a molecule can be altered, by transferring charge or applying an external field. Encouraging experimental indications of a conformationally mediated I-V have been obtained, for a fullerene-based transistor , and for the redox sidegroup-specific  or vibrationally mediated  NDR measurements. We are currently investigating the role of conformational changes, in conjunction with charging and gate electrostatics, in modulating the molecular I-V  in a three-terminal device.

(ii) Incoherence. For long molecules, molecular vibrations are important as sources of incoherent or inelastic scattering, leading to hopping or polaronic transport (see for e.g. ). Such inelastic effects can be naturally included in the NEGF formalism, requiring us to introduce another self-energy matrix describing the connection of the molecule with the source of the inelastic scattering (a phonon bath, for example) . In the tunneling regime, in particular, inelastic scattering turns out to be crucial for transport and dissipation in long molecular wires, such as DNA.

(iii) Correlation. Finally, there are examples of molecular measurements such as the Kondo effect , where the physics of current flow cannot be captured in terms of a simple one-particle picture, and require the incorporation of many-body correlation effects in our model. We leave these problems for future work.

We would like to thank P. Damle, D. Janes, A. Liang, S. Lodha, M. Lundstrom, M. Paulsson, T. Rakshit, R. Reifenberger, R. Venugopal and F. Zahid for useful discussions. This work has been supported by the US Army Research Office (ARO) under grant number DAAD19-99-1-0198.