text:
# Current rectification by simple molecular quantum dots: an ab-initio study

## I Introduction

Effective current rectifier is a necessary element of electronics circuitry. Aviram and Ratner suggested in 1974 that in a molecule containing donor ($`D)`$ and acceptor ($`A)`$ groups separated by a saturated $`\sigma `$-bridge (insulator) group, the (inelastic) electron transfer will be more favorable from $`A`$ to $`D,`$ rather than in the opposite direction. It was noted more recently that the electron transfer from $`D`$ to $`A`$ should start at smaller bias voltage because the highest occupied molecular orbital (HOMO) is centered on $`D`$ group, whereas the lowest unoccupied molecular orbital (LUMO) is located on the acceptor group $`A`$. The onset of resonant tunneling in this case corresponds to the alignment of HOMO on the $`D`$ group with LUMO on the $`A`$ group under external bias voltage, which are closest in energy. The molecular rectifiers (MR) synthesized in the laboratory were of somewhat different $`D\pi A`$ type, i.e. the “bridge” group was conjugated . The molecule, C16H$`{}_{33}{}^{}\gamma `$Q3CNQ, can be viewed as the fused naphthalene and TCNQ molecules with attached C16H32 alkane “tail” needed to provide enough van-der-Waals attraction between the molecules to assemble them into Langmuir-Blodgett film on a water. Although the molecule did show a rectification (with considerable hysteresis), it performed not like a conductor but rather like an anisotropic insulator because of that large alkane “tail”. Indeed, the currents reported in  were on the order of $`10^{17}`$ A/molecule, which make its application impractical. It was recently realized that in this molecule the resonance does not come from the alignment of HOMO and LUMO, since they cannot be decoupled through the conjugated $`\pi `$bridge, but rather due to an asymmetric voltage drop across the molecule where HOMO and LUMO are asymmetrically positioned with respect to the Fermi level of the electrodes . It is not clear if this mechanism has been observed. For instance, the contact between the alkane chain and electrode was certainly poor, via weak van-der-Waals forces. No temperature dependence of the resistivity was reported, but its large value may suggest that inelastic tunneling processes might have been involved.

It is clear that the rectifying molecular films built with the use of LB technique relying on use of long aliphatic (or other insulating) “tails” will produce hugely resistive molecules, not suitable for moletronics applications. There are reports of rectifying behavior in other classes of molecules, e.g. molecules chemisorbed on electrodes, where an observed asymmetry of I-V curves is likely due to asymmetric contacts with electrodes , or asymmetry of the molecule itself .

To achieve a good rectification by a molecule with a reasonable conductivity, one should avoid using molecules with long saturated (insulating) groups, which would make them prohibitively resistive. Thus, a suggestion to use relatively short molecules with certain end groups to allow their self-assembly on a metallic electrode’s surface seems to be very attractive. In that previous paper the transport through a phenyl ring connected to electrodes by asymmetric alkane chains (CH$`{}_{2}{}^{})_n`$ has been calculated with the use of a semi-empirical tight-binding method. Since the microscopic parameters of the model are poorly known, especially hopping integrals between molecule and electrodes, electrode work function and the affinity of the molecule (which may be strongly affected by bonding to a metal), the calculations have been performed for a number of these parameters. The calculations indicated a rectification ratio for -S-(CH$`{}_{2}{}^{})_2`$-C6H$`{}_{4}{}^{}`$(CH$`{}_{2}{}^{})_n`$S- molecule of about 100 for $`n=10`$ with the resistance $`R13`$G$`\mathrm{\Omega }.`$ It is difficult to know how quantitative these numbers are, because of a limited basis set and uncertainty in parameters of the model. For instance, an account for larger basis set for wave functions may lead to an increased estimate of the tunnel current, since the system would have more channels to conduct current. Experimentally, it may be easier to synthesize somewhat different compounds with a naphthalene group as a “molecular quantum dot”, like $``$S-(CH$`{}_{2}{}^{})_m`$-C10H$`{}_{6}{}^{}`$(CH$`{}_{2}{}^{})_n`$S$``$, which we will also refer to as $``$S-(CH$`{}_{2}{}^{})_m`$-Naph$``$(CH$`{}_{2}{}^{})_n`$S$`.`$ To obtain an accurate description of transport in this case, we employ an ab-initio non-equilibrium Green’s function method.

The calculation shows that indeed the current rectification $`I_+/I_{}100`$ may be possible for some designs like $``$S-(CH$`{}_{2}{}^{})_3`$-Naph$``$(CH$`{}_{2}{}^{})_{10}`$S$``$, where $`I_{+()}`$ is the forward (reverse) current. The difference between the forward and reverse voltages is, however, limited by other orbitals intervening into the conduction process. One needs the conducting orbital to be much closer in energy to the electrode Fermi level than the other ones (e.g. LUMO versus HOMO) and this energy asymmetry can be manipulated by “doping” the conjugated conducting part by attaching a donor (or acceptor) group, as will be shown below.

The present calculation takes into account only elastic tunneling processes. Inelastic processes may substantially modify results in the case of strong interaction of electrons with molecular vibrations, see . There are indications in the literature that the carrier might be trapped in a polaron state in saturated molecules slightly longer than we consider in the present paper.

## II Rectification by molecular quantum dot

The structure of the present molecular rectifier is shown in Fig. 1. The molecule consists of a central conjugated part (naphthalene) isolated from the electrodes by two insulating “arms” built from saturated aliphatic chains (CH2)n with lengths $`L_1`$ ($`L_2`$) for the left (right) chain. The principle of molecular rectification by a molecular quantum dot is illustrated in Fig. 2, where the electrically “active” molecular orbital, localized on the middle conjugated part, is LUMO, which lies at the energy $`\mathrm{\Delta }`$ above the electrode Fermi level at zero bias. The position of the LUMO, for instance, is determined by the work function of the metal $`q\varphi `$ and the affinity of the molecule $`q\chi ,`$ $`\mathrm{\Delta }=\mathrm{\Delta }_{LUMO}=q\left(\varphi \chi \right).`$ The position of HOMO is given by $`\mathrm{\Delta }_{HOMO}=\mathrm{\Delta }_{LUMO}E_g,`$ where $`E_g`$ is the HOMO-LUMO gap. If this orbital is considerably closer to the electrode Fermi level $`E_F,`$ then it will be brought into resonance with $`E_F`$ prior to other orbitals when the molecule is biased in either forward or reverse direction (Fig. 2). It is easy to estimate the forward and reverse bias voltage assuming that the voltage mainly falls off at the saturated (wide band gap) parts of the molecule with the lengths $`L_1`$ ($`L_2)`$ for the left (right) barrier, Fig. 2. Since the voltage drops proportionally to the lengths of the barriers, $`V_{1(2)}L_{1(2)},`$ we obtain for the partial voltage drops $`V_{1(2)}=VL_{1(2)}/L`$, where $`V`$ is the bias voltage. Forward and reverse voltages are found from the resonance condition, which gives

$`V_F`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{q}}(1+\xi ),`$ (1)
$`V_R`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{q}}\left(1+{\displaystyle \frac{1}{\xi }}\right),`$ (2)
$`V_F/V_R`$ $`=`$ $`\xi L_1/L_2,`$ (3)
where $`q`$ is the elementary charge. A significant difference between forward and reverse currents should be observed in the voltage range $`V_F<|V|<V_R`$. For our systems in question the parameter $`\mathrm{\Delta }=1.21.5`$ eV, Fig. 3,

the value determined by the work function of the electrodes and electron affinity of the conjugated part of the molecule, as well as by the strength of a surface dipole layer. For unlike electrodes with the work functions $`q\varphi _1`$ and $`q\varphi _2`$ there will be a built-in potential in the system $`q\varphi _b=q\left(\varphi _2\varphi _1\right)`$, Fig. 2(d). In this case the expressions for the forward and reverse bias voltage are modified as follows$`:`$

$`V_F`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{\Delta }}{q}}+\varphi _b\right)(1+\xi ),`$ (4)
$`V_R`$ $`=`$ $`\varphi _b+{\displaystyle \frac{\mathrm{\Delta }}{q}}\left(1+{\displaystyle \frac{1}{\xi }}\right).`$ (5)
Thus, the operating bias “window” in this case

$`V_RV_F`$ $`=`$ $`\left(V_RV_F\right)_0\xi \varphi _b`$ (6)
$`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{q}}\left({\displaystyle \frac{1}{\xi }}\xi \right)\xi \varphi _b.`$ (7)
Therefore, it slightly ($`\xi <1`$) increases when $`\varphi _b<0`$ ($`\varphi _1>\varphi _2`$, left work function is larger$`),`$ and shrinks down when $`\varphi _b>0`$ (right work function is larger). One would like to use as asymmetric a molecule as possible to increase the operating voltage, and in this case the difference of the work functions of the electrodes progressively becomes less important.

The main transport features of the present system can be qualitatively explained within a model of resonant tunneling through a localized state (molecular orbital) in the barrier. The current can be written as

$$I=\frac{2q}{h}𝑑ET(E)\left[f(E)f(E+qV)\right],$$
(8)
where the transmission $`T(E)`$at zero temperature is given by (for two spin directions)

$$T(E)=\frac{\mathrm{\Gamma }_L\mathrm{\Gamma }_R}{\left(EE_{MO}\right)^2+\left(\mathrm{\Gamma }_L+\mathrm{\Gamma }_R\right)^2/4},$$
(9)
where $`E_{MO}=E_F+\mathrm{\Delta }qV_1=E_F+\mathrm{\Delta }qVL_1/L`$ is the energy of the molecular orbital, relative to the Fermi level of the left electrode. The transmission is maximal and equals $`1`$ when $`E=E_{MO}`$ and $`\mathrm{\Gamma }_L=\mathrm{\Gamma }_R,`$ which corresponds to a symmetric position of the central conjugated part with respect to the electrodes. We can estimate $`\mathrm{\Gamma }_{L(R)}t^2/D=\mathrm{\Gamma }_0e^{2\kappa L_{1(2)}},`$ where $`t`$ is the overlap integral between the MO and the electrode, $`D`$ is the electron band width in the electrodes, $`\kappa `$ the inverse decay length of the resonant MO into the barrier. The latter quantity can be estimated from $`\mathrm{}^2\kappa ^2/2m^{}=U_0,`$ where $`U_0`$ is the barrier height \[$`4.8`$ eV in alkane chains (CH$`{}_{2}{}^{})_n`$,  $`5`$ eV from DFT calculations\], and $`m^{}1`$ the effective tunneling mass. The current above the threshold (when the resonant tunneling condition is met) is given by

$$I\frac{2q}{\mathrm{}}\mathrm{\Gamma }_0e^{2\kappa L_2}.$$
(10)
We see immediately that increasing the spatial asymmetry of the design (i.e. increasing $`L_2/L_1)`$ changes the operating voltage range linearly, $`V_RV_F\left(\mathrm{\Delta }/q\right)L_2/L_1,`$ and brings about an exponential decrease in current. This severely limits the ability to increase the rectification ratio while simultaneously keeping the resistance of the molecule at a reasonable value. One therefore needs a quantitative analysis of all microscopic details, including molecule-electrode contact, and potential and charge distribution, in order to make comparisons with the data.

## III Method

Since we study relatively short molecules, electron transport is likely to proceed elastically. We shall treat the elastic tunneling processes from first principles in order to evaluate all the relevant microscopic characteristics within the density functional theory. This allows us to check the accuracy of semi-empirical methods used previously, and find results for new objects of experimental interest, like the present molecular quantum dot with a naphthalene conjugated middle group. We use an ab-initio approach that combines the Keldysh non-equilibrium Green’s function (NEGF)  with pseudopotential-based real space density functional theory (DFT). In contrast to standard DFT work, our NEGF-DFT approach considers a quantum mechanical system with an open boundary condition provided by semi-infinite electrodes under external bias voltage. The main advantages of our approach are (i) a proper treatment of the open boundary condition; (ii) a fully atomistic treatment of the electrodes and (iii) a self-consistent calculation of the non-equilibrium charge density using NEGF.

For the open devices studied here, our system is divided into three regions, a left and right electrode and a central scattering region, which contains a portion of our semi-infinite electrodes (see Fig. 1). By including atomic layers of the electrode in our scattering region, we provide a proper treatment of the electrode-molecule interaction, and allow for the effective screening of the molecule away from our molecular junction. External bias voltage can then be applied as a boundary condition for the Poisson equation in our central scattering region. The Hartree potential in the scattering region is solved using a multigrid real-space numerical approach. Finally, the non-equilibrium density matrix $`\widehat{\rho }`$ is constructed via NEGF ,

$$\widehat{\rho }=\frac{1}{2\pi i}𝑑E\mathrm{Tr}G^<(E),$$
(11)
where the lesser Green’s function is determined from the operator relation

$$G^<=G^R\mathrm{\Sigma }^<G^A,$$
(12)
via the retarded (advanced) Green’s functions $`G^{R(A)}`$ and $`\mathrm{\Sigma }^<`$ the self-energy part describing an injection of charge into the scattering region (molecule) from the electrodes. This process is characterized by the usual open channel representation of the scattering states with momenta $`k_l^n`$ ($`k_r^n)`$ in the left (right) electrodes, where index $`n`$ enumerates all open channels (available Bloch states of the electrons in the electrodes). The Fermi functions $`f(k_l^n)`$ and $`f(k_r^n)`$ define the occupied states in the leads to be accounted for in evaluation of the electron transport. To evaluate the Green’s function, one treats the disconnected system consisting of the left (right) lead marked by index $`l`$ $`(r)`$ and the scatterer $`c.`$ The transport Green’s function is then found from the Dyson equation

$$\left(G^R\right)^1=\left(G_0^R\right)^1V,$$
(13)
where the unperturbed retarded Green’s function is defined in operator form as

$`\left(G_0^R\right)^1`$ $`=`$ $`(E+i0)\widehat{S}\widehat{H},`$ (14)
$`\widehat{H}`$ $`=`$ $`\left(\begin{array}{ccc}H_{l,l}& H_{l,c}& 0\\ H_{c,l}& H_{c,c}& H_{c,r}\\ 0& H_{r,c}& H_{r,r}\end{array}\right).`$ (18)
$`H`$ is the Hamiltonian matrix, $`H_{c,l}=H_{l,c}^{}`$, $`H_{r,c}=H_{c,r}^{}`$, $`S`$ is the overlap matrix, $`S_{i,j}=\chi _i|\chi _j`$ for non-orthogonal basis set orbitals $`\chi _i`$, and the coupling of the scatterer to the leads is given by the Hamiltonian matrix $`V=\mathrm{diag}[\mathrm{\Sigma }_{l,l},0,\mathrm{\Sigma }_{r,r}].`$ In all formulas above it is assumed that matrix elements are given in terms of the non-orthogonal basis set.

The self-energy part $`\mathrm{\Sigma }^<,`$ which is used to construct the non-equilibrium electron density in the scattering region, is found from

$$\mathrm{\Sigma }^<=2iIm\left[f(E)\mathrm{\Sigma }_{l,l}+f(E+qV)\mathrm{\Sigma }_{r,r}\right].$$
(19)
This expression accounts for the steady charge “flowing in” from the electrodes. The transmission probability, which determines the current according to Eq. (8), is given by

$$T(E)=4\mathrm{T}\mathrm{r}\left[Im\left(\mathrm{\Sigma }_{l,l}\right)G_{l,r}^RIm\left(\mathrm{\Sigma }_{r,r}\right)G_{r,l}^A\right].$$
(20)
The local density of states, which is a very important quantity to characterize the spatial character of the wave functions, can be obtained from

$$N(E,\stackrel{}{r})=\frac{1}{2\pi i}\mathrm{Tr}\left[G^<(E)|\chi (\stackrel{}{r})\chi (\stackrel{}{r})|\right].$$
(21)
To construct the system Hamiltonian, we define the atomic cores using standard norm-conserving nonlocal pseudopotentials, and expand the Kohn-Sham wave functions in an s-,p-,d- real space atomic orbital basis . External bias is incorporated into the Hartree potential, and in this way the non-equilibrium charge density is iterated to self-consistency, after which the current-voltage (I-V) characteristics have been computed.

## IV Results

The physical picture of the transmission through a series of rectifiers $``$S-(CH$`{}_{2}{}^{})_m`$-C10H$`{}_{6}{}^{}`$(CH$`{}_{2}{}^{})_n`$S$``$ for $`m=2`$ and $`n=2,4,6,10`$ is illustrated in Figs. 3,4. The transmission coefficient $`T(E)`$ (top panel) and the density of states $`N(E)=𝑑VN(E,\stackrel{}{r})`$ (bottom panel) are shown at zero bias voltage in Fig. 3. The Fermi level in the electrodes is taken as the energy origin. It follows from this figure that indeed the LUMO is the closest molecular orbital transparent to electron transport, and it is higher in energy than $`E_F`$ by an amount $`\mathrm{\Delta }=1.21.5`$ eV depending on the molecule. Since the molecule is asymmetric, even at the resonance energy we have $`T(E_{LUMO})0.5,`$ Fig. 3.

This can be traced back to the Breit-Wigner formula (9), which for $`\mathrm{\Gamma }_R<\mathrm{\Gamma }_L`$ gives $`T(E_{LUMO})4\mathrm{\Gamma }_R/\mathrm{\Gamma }_L<1,`$ and it falls off as the asymmetry increases. There the transmission through nominal HOMO state is negligible, but HOMO-2 state conducts very well. As we shall see below, the HOMO-2 defines the threshold reverse voltage $`V_R,`$ thus limiting the operating voltage range (7).

The maps of the densities of states $`N(E_{LUMO},\stackrel{}{r})`$ in Fig. 4 show the development of the tunneling gap between the naphthalene group and more distant electrode separated by alkane chain (CH2)$`_n,`$ $`n=210.`$ The symmetry and weight distribution of LUMO in the device is very similar to LUMO in the isolated molecule, telling us that the character of the molecular wave functions is preserved to a high degree. Obviously, in contact with electrodes molecular orbitals acquire some finite width $`\mathrm{\Gamma }=\mathrm{\Gamma }_L+\mathrm{\Gamma }_R,`$ and for highly asymmetric molecules $`\mathrm{\Gamma }\mathrm{\Gamma }_L`$, so it is mainly defined by the transparency of the shorter (left) insulating group (CH2)m. One needs a rather small width $`\mathrm{\Gamma }`$ which defines the sharpness of the voltage threshold for the current turn-on . We have calculated the molecules with $`m=0,1,2`$ with the results for $`\mathrm{\Gamma }_m=130,`$ $`60`$, and 12 meV (Fig. 5). The last value for the left insulating group (CH$`{}_{2}{}^{})_2`$ seems to be most reasonable, since the resonant peak is rather narrow yet the molecule remains conductive. However, in some cases the short insulating group (CH$`{}_{2}{}^{})_3`$ may be preferential. This conclusion is similar to the one reached on the basis of the tight-binding model, but the peak positions and size are substantially different in both cases, the peak size being much larger in the case of the present DFT calculations, Fig. 5.

Our assumption that the voltage drop is proportional to the lengths of the alkane groups on both sides, is quantified by the calculated potential ramp in Fig. 6.

It is close to a linear slope along the (CH$`{}_{2}{}^{})_n`$ chains. There is also a smaller but still noticeable voltage drop across the naphthalene group, which can be viewed as a fragment of wide-band semiconductor.

The threshold bias voltages can be obtained from Fig. 7,

showing the evolution of the LUMO versus bias voltage with respect to the chemical potentials of the right (left) electrodes $`\mu _{R(L)}=E_F+()qV/2.`$ The forward voltage corresponds to the crossing of the LUMO($`V)`$ and $`\mu _R(V)`$, which happens at about 2V. Extrapolated crossing of LUMO and $`\mu _L`$ is at large negative voltages but, unfortunately, this large desirable difference between $`V_F`$ and $`V_R`$ does not materialize. Although the LUMO defines the forward threshold voltages in all the molecules studied here, the reverse voltage is defined by the HOMO-2 for “right” barriers (CH$`{}_{2}{}^{})_n`$ with $`n=6,10`$. The I-V curves are plotted in Fig. 8

and the corresponding parameters are listed in the Table. We see that the rectification ratio for current in the operation window $`I_+/I_{}`$ reaches the maximum value of 35 for “2-10” molecule ($`m=2,`$ $`n=10).`$ It was estimated to be about 100 for “2-10” molecule with a phenyl ring as a conjugated central part. We have checked a series of molecules with central phenyl ring and, unfortunately, do not find any significant rectification in this case, because the Fermi level in DFT calculations appears to lie close to mid-gap. Thus, the tight-binding results for large rectification in one-phenyl ring molecules might have been the results of the model. The estimated resistance of the molecule in DFT is substantially less (2.8 G$`\mathrm{\Omega })`$ compared to the tight-binding model with a phenyl ring as a conjugated central unit (13 G$`\mathrm{\Omega })`$, see Table.

One can manipulate the system in order to increase the energy asymmetry of conducting orbitals (reduce $`\mathrm{\Delta })`$. To shift the LUMO towards $`E_F`$, one needs in effect to gate the naphthalene group by positive voltage. A chemical way of doing this would be to attach some chemical group, like -C$``$N in the present case. This group withdraws a fraction of an electron and acquires the negative charge $`\delta =`$ $`0.11q.`$ The charge is donated mainly by the naphthalene group, which acquires a positive charge and orbitals on the molecule shift to lower energies almost rigidly by about 0.6 eV, Fig. 9.

On the resulting I-V curve, both thresholds $`V_{F(R)}`$ now correspond to LUMO crossing either the right or left electrode Fermi levels (see Fig. 7), and the operating voltage range is modified significantly.

## V Conclusions

We have presented parameter-free DFT calculations of a class of molecular quantum dots showing a current rectification at a reasonable current density through the molecule. The explored mechanism of rectification is different from Aviram-Ratner diodes and relies on considerable asymmetry of the spatial composition and energy structure of conducting orbitals of the molecule. Since there is no empirical parameters, the predictions are quantitative and can be tested in a controllable fashion. We have studied the molecular quantum dots with the naphthalene central conjugated group, and found a rectification of about 35 for $``$S-(CH$`{}_{2}{}^{})_2`$-C10H$`{}_{6}{}^{}`$(CH$`{}_{2}{}^{})_6`$S$``$ molecules. The modest rectification ratio is a consequence of the resonant processes at these molecular sizes. The present study suggests also that the rectification ratio in the case of a central single phenyl ring, studied previously, is dramatically overestimated by the tight-binding model.

The rectification ratio is not great by any means, but one should bear in mind that this is a device necessarily operating in a ballistic quantum-mechanical regime, because of small size. This is very different from present Si devices with carriers diffusing through the system. As they become smaller, however, the same effects, as those discussed here for molecular rectifiers, will eventually take over and will tend to diminish the rectification ratio. In this regard the present results are indicative of the problems which will ultimately be faced by Si devices of molecular proportions.

It is important that the parameters of the device, like the threshold voltages for forward and reverse bias, can be significantly modified by “chemical” doping. In the present case adding the -C$``$N electron withdrawing group has scaled the threshold range by a factor of about two. There may be other modifications that could improve/add functionality of the molecular devices by varying central, end, and side groups of the molecule. The present results illustrate very clearly that the optimization of the parameters of the molecular devices is typically facing the problem that some parameters, like level widths, vary exponentially, while others, like the voltage thresholds, vary as a power law– as a function of the molecular design. Although the approximate band diagrams illustrating the level structure, Fig. 2, are very helpful for qualitative description, the atomistic details of the electron density distribution are extremely important, Fig. 4. This necessitates the use of ab-initio methods, without which one may obtain qualitatively wrong results.

We have found that in the present molecules significant rectification may be obtained in relatively short self-assembled molecules. Self-assembly provides for good contact with electrode and relatively small resistance per molecule. In existing studies of Aviram-Ratner-type molecules they were assembled by Langmuir-Blodgett deposition and are prohibitively resistive by design because of a long aliphatic tail needed for the assembly . Thus, for applications the self-assembly of shorter molecules on metallic electrodes seems to be essential. One should keep in mind that the molecular rectification, here estimated for ideal contacts at zero temperature, will be reduced at finite temperature and due to inevitable disorder in the contact area and in the molecular film. Indeed, the optimal width of the level was found to be about 12 meV, which is just a half of the room temperature, so the temperature (and static disorder) will significantly broaden the level and reduce the rectification. In terms of simplest cross-point memories one may use a large number of molecules on a pitch. This may allow for use of longer, better rectifying molecules. There, however, the inelastic processes will play a progressively more significant role and will further reduce the rectification. Further studies of molecular rectifiers should be focused on these outstanding issues in conjunction with experimental studies.

We acknowledge useful discussions with Shun-chi Chang and R. S. Williams. The work has been supported by DARPA.