text:
# Resistance of a Molecule

## 1 Introduction

In recent years, several experimental groups have reported measurements of the current-voltage (I-V) characteristics of individual or small numbers of molecules. Even three-terminal measurements showing evidence of transistor action has been reported using carbon nanotubes  as well as self-assembled monolayers of conjugated polymers . These developments have attracted much attention from the semiconductor industry and there is great interest from an applied point of view to model and understand the capabilities of molecular conductors. At the same time, this is also a topic of great interest from the point of view of basic physics. A molecule represents a quantum dot, at least an order of magnitude smaller than semiconductor quantum dots, which allows us to study many of the same mesoscopic and/or many-body effects at far higher temperatures.

So what is the resistance of a molecule? More specifically, what do we see when we connect a short molecule between two metallic contacts as shown in Fig. 1 and measure the current (I) as a function of the voltage (V)? Most commonly we get I-V characteristics of the type sketched in Fig. 2. This has been observed using many different approaches including breakjunctions , scanning probes , nanopores  and a host of other methods (see for example ). A number of theoretical models have been developed for calculating the I-V characteristics of molecular wires using semi-empirical  as well as first principles  theory.

Our purpose in this chapter is to provide an intuitive explanation for the observed I-V characteristics using simple models to illustrate the basic physics. However, it should be noted that molecular electronics is a fast developing field and much of the excitement arises from the possibility of discovering novel physics beyond the paradigms discussed here. To cite a simple example, very few experiments to date  incorporate the gate electrode shown in Fig. 1, and we will largely ignore the gate in this chapter. However, the gate electrode can play a significant role in shaping the I-V characteristics and deserves more attention. This is easily appreciated by looking at the applied potential profile $`U_{app}`$ generated by the electrodes in the absence of the molecule. This potential profile satisfies the Laplace equation without any net charge anywhere, and is obtained by solving:

$$\left(ϵU_{app}\right)=0$$
(1)
subject to the appropriate boundary values on the electrodes (Fig. 3). It is apparent that the electrode geometry has a significant influence on the potential profile that it imposes on the molecular species and this in term could obviously affect the I-V characteristics in a significant way.

After all, it is well known that a three-terminal metal/oxide/semiconductor Field Effect Transistor (MOSFET) with a gate electrode has a very different I-V characteristic compared to a two terminal “n-i-n” diode: The current in a MOSFET saturates under increasing bias, but the current in an “n-i-n” diode keeps increasing indefinitely. In contrast to the MOSFET, whose I-V is largely dominated by classical electrostatics, the I-V characteristics of molecules is determined by a more interesting interplay between nineteenth century physics (electrostatics) and twentieth century physics (quantum transport) and it is important to do justice to both aspects.

*Outline of chapter*: We will start in section 2 with a qualitative discussion of the main factors affecting the I-V characteristics of molecular conductors, using a simple toy model to illustrate their role. However, this toy model misses two important factors: (1) Shift in the energy level due to *charging effects* as the molecule loses or gains electrons and (2) *broadening* of the energy levels due to their finite lifetime arising from the coupling ($`\mathrm{\Gamma }_1`$ and $`\mathrm{\Gamma }_2`$) to the two contacts. Once we incorporate these effects (section 3) we obtain more realistic I-V plots, even though the toy model assumes that conduction takes place independently through individual molecular levels. In general, however, multiple energy levels are simultaneously involved in the conduction process. In section 4 we will describe the non-equilibrium Green’s function (NEGF) formalism which can be viewed as a generalized version of the one-level model to include multiple levels or conduction channels. This formalism provides a convenient framework for describing quantum transport  and can be used in conjunction with ab initio or semi-empirical Hamiltonians as described in a set of related articles . Here (section 5) we will illustrate the NEGF formalism with a simple semi-empirical model for a gold wire, ’$`n`$’ atoms long and one atom in cross-section. We could call this a $`\text{Au}_n`$ molecule through that is not how one normally thinks of a gold wire. However, this example is particularly instructive because it shows the lowest possible “resistance of a molecule” per channel, which is $`\pi \mathrm{}/e^2=12.9\text{k}\mathrm{\Omega }`$ .

## 2 Qualitative Discussion

### 2.1 Where is the Fermi energy?

*Energy Level Diagram*: The first step in understanding the current (I) vs. voltage (V) curve for a molecular conductor is to draw an energy level diagram and locate the Fermi energy. Consider first a molecule sandwiched between two metallic contacts, but with very weak electronic coupling. We could then line up the energy levels as shown in Fig. 4 using the metallic work function (WF) and the electronic affinity ($`EA`$) and ionization potential ($`IP`$) of the molecule. For example, a (111) gold surface has a work function of $`5.3`$ eV while the electron affinity and ionization potential, $`EA_0`$ and $`IP_0`$, for isolated phenyl dithiol (Fig. 1) in the gas phase have been reported to be $`2.4`$ eV and $`8.3`$ eV respectively . These values are associated with electron emission and injection to and from a vacuum and may need some modification to account for the metallic contacts. For example the actual $`EA`$, $`IP`$ will possibly be modified from $`EA_0`$, $`IP_0`$ due to the image potential $`W_{im}`$ associated with the metallic contacts :

$`EA`$ $`=`$ $`EA_0+W_{im}`$ (2)
$`IP`$ $`=`$ $`IP_0W_{im}`$ (3)
The probability of the molecule losing an electron to form a positive ion is equal to $`e^{(WFIP)/k_BT}`$ while the probability of the molecule gaining an electron to form a negative ion is equal to $`e^{(EAWF)/k_BT}`$. We thus expect the molecule to remain neutral as long as both $`(IPWF)`$ and $`(WFEA)`$ are much larger than $`k_BT`$, a condition that is usually satisfied for most metal-molecule combinations. Since it costs too much energy to transfer one electron into or out of the molecule, it prefers to remain neutral in equilibrium.

The picture changes qualitatively if the molecule is chemisorbed directly on the metallic contact (Fig. 5). The molecular energy levels are now broadened significantly by the strong hybridization with the delocalized metallic wavefunctions, making it possible to transfer fractional amounts of charge to or from the molecule. Indeed there is a charge transfer which causes a change in the electrostatic potential inside the molecule and the energy levels of the molecule are shifted by a contact potential (CP), as shown.

It is now more appropriate to describe transport in terms of the HOMO-LUMO levels associated with incremental charge transfer  rather than the affinity and ionization levels associated with integer charge transfer. Whether the molecule-metal coupling is strong enough for this to occur depends on the relative magnitudes of the single electron charging energy ($`U`$) and energy level broadening ($`\mathrm{\Gamma }`$). As a rule of thumb, if $`U>>\mathrm{\Gamma }`$, we can expect the structure to be in the Coulomb Blockade (CB) regime characterized by integer charge transfer; otherwise it is in the self-consistent field (SCF) regime characterized by fractional charge transfer. This is basically the same criterion that one uses for the Mott transition in periodic structures, with $`\mathrm{\Gamma }`$ playing the role of the hopping matrix element. It is important to note that for a structure to be in the CB regime both contacts must be weakly coupled, since the total broadening $`\mathrm{\Gamma }`$ is the sum of the individual broadening due to the two contacts. Even if only one of the contacts is coupled strongly we can expect $`\mathrm{\Gamma }U`$ thus putting the structure in the SCF regime. Fig. 14 in Sec. 3.2 illustrates the I-V characteristics in the CB regime using a toy model. However, a moderate amount of broadening destroys this effect (see Fig. 17) and in this chapter we will generally assume that the conduction is in the SCF regime.

*Location of the Fermi energy:* The location of the Fermi energy relative to the HOMO and LUMO levels is probably the most important factor in determining the current (I) versus voltage (V) characteristics of molecular conductors. Usually it lies somewhere inside the HOMO-LUMO gap. To see this, we first note that $`E_f`$ is located by the requirement that the number of states below the Fermi energy must be equal to the number of electrons in the molecule. But this number need not be equal to the integer number we expect for a neutral molecule. A molecule does not remain exactly neutral when connected to the contacts. It can and does pick up a fractional charge depending on the work function of the metal. However, the charge transferred ($`\delta n`$) for most metal-molecule combinations is usually much less than one. If $`\delta n`$ were equal to +1, the Fermi energy would lie on the LUMO while if $`\delta n`$ were -1, it would lie on the HOMO. Clearly for values in between, it should lie somewhere in the HOMO-LUMO gap.

A number of authors have performed detailed calculations to locate the Fermi energy with respect to the molecular levels for a phenyl dithiol molecule sandwiched between gold contacts, but there is considerable disagreement. Different theoretical groups have placed it close to the LUMO  or to the HOMO . The density of states inside the HOMO-LUMO gap is quite small making the precise location of the Fermi energy very sensitive to small amounts of electron transfer, a fact that could have a significant effect on both theory and experiment. As such it seems justifiable to treat $`E_f`$ as a “fitting parameter” within reasonable limits when trying to explain experimental I-V curves.

*Broadening by the contacts*: “Common sense” suggests that the strength of coupling of the molecule to the contacts is important in determining the current flow - the stronger the coupling, the larger the current. A useful quantitative measure of the coupling is the resulting broadening $`\mathrm{\Gamma }`$ of the molecular energy levels, see Fig. 6. This broadening $`\mathrm{\Gamma }`$ can also be related to the time $`\tau `$ it takes for an electron placed in that level to escape into the contact: $`\mathrm{\Gamma }=\mathrm{}/\tau `$. In general, the broadening $`\mathrm{\Gamma }`$ could be different for different energy levels. Also it is convenient to define two quantities $`\mathrm{\Gamma }_1`$ and $`\mathrm{\Gamma }_2`$, one for each contact, with the total broadening $`\mathrm{\Gamma }=\mathrm{\Gamma }_1+\mathrm{\Gamma }_2`$.

One subtle point. Suppose an energy level is located well below the Fermi energy in the contact, so that the electrons are prevented from escaping by the exclusion principle. Would $`\mathrm{\Gamma }`$ be zero? No, the broadening would still be $`\mathrm{\Gamma }`$, independent of the degree of filling of the contact as discussed in Ref. . This observation is implicit in the NEGF formalism, though we do not invoke it explicitly.

### 2.2 Current flow as a “balancing act”

Once we have drawn an equilibrium energy level diagram, we can understand the process of current flow which involves a non-equilibrium situation where the different reservoirs (e.g., the source and the drain) have different electrochemical potentials $`\mu `$, see Fig. 7. For example if a positive voltage V is applied externally to the drain with respect to the source, then the drain has an electrochemical potential lower than that of the source by $`eV`$: $`\mu _2=\mu _1eV`$. The source and drain contacts thus have different Fermi functions and each seeks to bring the active device into equilibrium with itself. The source keeps pumping electrons into it hoping to establish equilibrium. But equilibrium is never achieved as the drain keeps pulling electrons out in its bid to establish equilibrium with itself. The device is thus forced into a balancing act between two reservoirs with different agendas which sends it into a non-equilibrium state intermediate between what the source and drain would like to see. To describe this balancing process we need a kinetic equation that keeps track of the in and out-flow of electrons from each of the reservoirs.

*Kinetic equation*: This balancing act is easy to see if we consider a simple one level system, biased such that the energy $`ϵ`$ lies in between the electrochemical potentials of the two contacts (Fig. 8). An electron in this level can escape into contacts 1 and 2 at a rate of $`\mathrm{\Gamma }_1/\mathrm{}`$ and $`\mathrm{\Gamma }_2/\mathrm{}`$ respectively. If the level were in equilibrium with contact 1 then the number of electrons occupying the level would be given by:

$$N_1=2(\text{for spin})f(ϵ,\mu _1)$$
(4)
where

$$f(ϵ,\mu )=\frac{1}{1+e^{\frac{ϵ\mu }{k_BT}}}$$
(5)
is the Fermi function. Similarly if the level were in equilibrium with contact 2 the number would be:

$$N_2=2(\text{for spin})f(ϵ,\mu _2)$$
(6)
Under non-equilibrium conditions the number of electrons N will be somewhere in between $`N_1`$ and $`N_2`$. To determine this number we write a steady state *kinetic equation* that equates the net current at the left junction:

$$I_L=\frac{e\mathrm{\Gamma }_1}{\mathrm{}}(N_1N)$$
(7)
to the net current at the right junction:

$$I_R=\frac{e\mathrm{\Gamma }_2}{\mathrm{}}(NN_2)$$
(8)
Steady state requires $`I_L=I_R`$, from which we obtain:

$$N=2\frac{\mathrm{\Gamma }_1f(ϵ,\mu _1)+\mathrm{\Gamma }_2f(ϵ,\mu _2)}{\mathrm{\Gamma }_1+\mathrm{\Gamma }_2}$$
(9)
so that from Eq. 7 or 8 we obtain the current:

$$I=\frac{2e}{\mathrm{}}\frac{\mathrm{\Gamma }_1\mathrm{\Gamma }_2}{\mathrm{\Gamma }_1+\mathrm{\Gamma }_2}\left(f(ϵ,\mu _1)f(ϵ,\mu _2)\right)$$
(10)
Eq. 10 follows very simply from an elementary model, but it serves to illustrate a basic fact about the process of current flow. No current will flow if $`f(ϵ,\mu _1)=f(ϵ,\mu _2)`$. A level that is way below both electrochemical potentials $`\mu _1`$ and $`\mu _2`$ will have $`f(ϵ,\mu _1)=f(ϵ,\mu _2)=1`$ and will not contribute to the current, just like a level that is way above both potentials $`\mu _1`$ and $`\mu _2`$ and has $`f(ϵ,\mu _1)=f(ϵ,\mu _2)=0`$. It is only when the level lies in between $`\mu _1`$ and $`\mu _2`$ (or within a few $`k_BT`$ of $`\mu _1`$ and $`\mu _2`$) that we have $`f(ϵ,\mu _1)f(ϵ,\mu _2)`$ and a current flows. Current flow is thus the result of the *difference in opinion* between the contacts. One contact would like to see more electrons (than $`N`$) occupy the level and keeps pumping them in, while the other would like to see fewer than $`N`$ electrons and keeps pulling them out. The net effect is a continuous transfer of electrons from one contact to another.

Fig. 9 shows a typical current ($`I`$) versus voltage ($`V`$) calculated from Eq. 10, using the parameters indicated in the caption. At first the current is zero because both $`\mu _1`$, $`\mu _2`$ are above the energy level.

Once $`\mu _2`$ drops below the energy level, the current increases to $`I_{max}`$ which is the maximum current that can flow through one level and is obtained from Eq. 10 by setting $`f(ϵ,\mu _1)=1`$ and $`f(ϵ,\mu _2)=0`$:

$$I_{max}=\frac{2e}{\mathrm{}}\mathrm{\Gamma }_{\text{eff}}=\frac{2e}{\mathrm{}}\frac{\mathrm{\Gamma }_1\mathrm{\Gamma }_2}{\mathrm{\Gamma }_1+\mathrm{\Gamma }_2}$$
(11)
Note that in Fig. 9 we have set $`\mu _1=E_feV/2`$ and $`\mu _2=E_f+eV/2`$. We could, of coarse, just as well have set $`\mu _1=E_feV`$ and $`\mu _2=E_f`$. But the, the average potential in the molecule would be $`V/2`$ and we would need to shift $`ϵ`$ appropriately. It is more convenient to choose our reference such that the average molecular potential is zero and there is no need to shift $`ϵ`$.

Note that the current is proportional to $`\mathrm{\Gamma }_{\text{eff}}`$ which is the parallel combination of $`\mathrm{\Gamma }_1`$ and $`\mathrm{\Gamma }_2`$. This seems quite reasonable if we recognize that $`\mathrm{\Gamma }_1`$ and $`\mathrm{\Gamma }_2`$ represents the strength of the coupling to the two contacts and as such are like two *conductances in series*. For long conductors we would expect a third conductance in series representing the actual conductor. This is what we usually have in mind when we speak of conductance. But short conductors have virtually zero resistance and what we measure is essentially the contact or interface resistance111 Four-terminal measurments have been used to separate the contact from the device resistance (see for example Ref. ).. This is an important conceptual issue, that caused much argument and controversy in the 1980’s and was finally resolved when experimentalists measured the conductance of very short conductors and found it approximately equal to $`2e^2/h`$ which is a fundamental constant equal to 77.8 $`\mu `$A/V. The inverse of this conductance $`h/2e^2=12.9`$ k$`\mathrm{\Omega }`$ is now believed to represent the minimum contact resistance that can be achieved for a one-channel conductor. Even a copper wire with a one atom cross section will have a resistance at least this large. Our simple one-level model (Fig, 9) does not predict this result because we have treated the level as discrete, but the more complete treatment in later sections will show it.

## 3 Coulomb Blockade?

As we mentioned in section 2, a basic question we need to answer is whether the process of conduction through the molecule belongs to the Coulomb Blockade (CB) or the Self-Consistent Field (SCF) regime. In this section, we will first discuss a simple model for charging effects (section 3.1) and then look at the distinction between the simple SCF regime and the CB regime (section 3.2). Finally, in section 3.3 we show how moderate amount of level broadening often destroy CB effects making a simple SCF treatment quite accurate.

### 3.1 Charging Effects

Given the level ($`ϵ`$), broadening ($`\mathrm{\Gamma }_1`$, $`\mathrm{\Gamma }_2`$) and the electrochemical potentials $`\mu _1`$ and $`\mu _2`$ of the two contacts, we can solve Eq. 10 for the current $`I`$. But we want to include charging effects in the calculations. Therefore, we add a potential $`U_{SCF}`$ due to the change in the number of electrons from the equilibrium value ($`f_0=f(ϵ_0,E_f)`$):

$$U_{SCF}=U\left(N2f_0\right)$$
(12)
similar to a Hubbard model. We then let the level $`ϵ`$ float up or down by this potential:

$$ϵ=ϵ_0+U_{SCF}$$
(13)
Since the potential depends on the number of electrons, we need to calculate the potential using the self consistent procedure shown in Fig. 10.

Once the converged solution is obtained, the current is calculated from Eq. 10. This very simple model captures much of the observed physics of molecular conduction. For example, the results obtained by setting $`E_f=5.0`$ eV, $`ϵ_0=5.5`$ eV, $`\mathrm{\Gamma }_1=0.2`$ eV, $`\mathrm{\Gamma }_2=0.2`$ eV are shown in Fig. 11 with ($`U=1.0`$ eV) and without ($`U=0`$ eV) charging effects. The finite width of the conductance peak (with $`U=0`$) is due to the temperature used in the calculations ($`k_BT=0.025`$ eV). Note how the inclusion of charging tends to broaden the sharp peaks in conductance, even though we have not included any extra level broadening in this calculation. The size of the conductance gap is directly related to the energy difference between the molecular energy level and the Fermi energy. The current starts to increase when the voltage reaches $`1`$ V, which is exactly $`2\left|E_fϵ_0\right|`$ as would be expected even from a theory with no charging. Charging enters the picture only at higher voltages, when a chemical potential tries to cross the level. The energy level shifts in energy (Eq. 13) if the charging energy is non-zero. Thus, for a small charging energy, the chemical potential easily crosses the level giving a sharp increase of the current. If the charging energy is large, the current increases gradually since the energy level follows the chemical potential due to the charging.

What determines the conductance gap? The above discussion shows that the conductance gap is equal to $`4\left(\left|E_fϵ_0\right|\mathrm{\Delta }\right)`$ where $`\mathrm{\Delta }`$ is equal to $`4k_BT`$ (plus $`\mathrm{\Gamma }_1+\mathrm{\Gamma }_2`$ if broadening is included, see section 3.3), and $`ϵ_0`$ is the HOMO or LUMO level whichever is closest to the Fermi energy, as pointed out in Ref. . This is unappreciated by many who associate the conductance gap with the HOMO-LUMO gap. However, we believe that what conductance measurements show is the gap between the Fermi energy and the nearest molecular level222 With very asymmetric contacs, the conductance gap could be equal to the HOMO-LUMO gap as commonly assumed in interpreting STM spectra. However, we belive that the picture presented here is more accurate unless the contact is so strongly coupled that there is a significant density of Metal-Induced Gap States (MIGS). Fig. 12 shows the I-V characteristics calculated using a two-level model (obtained by a straightforward extension of the one-level model) with the Fermi energy located differently within the HOMO-LUMO gap giving different conductance gaps corresponding to the different values of $`\left|E_fϵ_0\right|`$. Note that with the Fermi energy located halfway in between, the conductance gap is twice the HOMO-LUMO gap and the I-V shows no evidence of charging effects because the depletion of the HOMO is neutralized by the charging of the LUMO. This perfect compensation is unlikely in practice, since the two levels will not couple identically to the contacts as assumed in the model.

A very interesting effect that can be observed is the asymmetry of the I-V characteristics if $`\mathrm{\Gamma }_1\mathrm{\Gamma }_2`$ as shown in Fig. 13. This may explain several experimental results which show asymmetric I-V  as discussed by Ghosh et al. . Assuming that the current is conducted through the HOMO level ($`E_f>ϵ_0`$), the current is less when a positive voltage is applied to the strongly coupled contact, see Fig. 13(a). This is due to the effects of charging as has been discussed in more detail in Ref. . Ghosh et al. also shows that this result will reverse if the conduction is through the LUMO level. We can simulate this situation by setting $`ϵ_0`$ equal to $`4.5`$ eV, $`0.5`$ eV above the equilibrium Fermi energy $`E_f`$. The sense of asymmetry is now reversed as shown in Fig. 13(b). The current is larger when a positive voltage is applied to the strongly coupled contact. Comparing with STM measurements seems to favor the first case, i.e., conduction through the HOMO .

### 3.2 Unrestricted Model

In the previous examples (Figs. 11, 13) we have used values of $`\mathrm{\Gamma }_{1,2}`$ that are smaller than the charging energy $`U`$. However, under these conditions one can expect Coulomb Blockade (CB) effects which are not captured by a “restricted solution” which assumes that both spin orbitals see the same self-consistent field. However, an unrestricted solution, which allows the spin degeneracy to be lifted, will show these effects333 The unrestricted one-particle picture discussed here provides at least a reasonable qualitative picture of CB effects, though a complete description requires a more advanced many particle picture . The one-particle picture leads to one of many possible states of the device depending on our initial guess, while a full many particle picture would include all states.. For example, if we replace Eq. 13 with ($`f_0=f(ϵ_0,E_f)`$):

$`ϵ_{}=ϵ_0+U\left(N_{}f_0\right)`$ (14)
$`ϵ_{}=ϵ_0+U\left(N_{}f_0\right)`$ (15)
where the up-spin level feels a potential due to the down-spin electrons and vice-versa, then we obtain I-V curves as shown in Fig. 14.

If the SCF iteration is started with a spin degenerate solution, the same restricted solution as before is obtained. However, if the iteration is started with a spin non-degenerate solution a different looking I-V is obtained. The electrons only interact with the the electron of the opposite spin. Therefore, the chemical potential of one contact can cross one energy level of the molecule since the charging of that level only affects the opposite spin level. Thus, the I-V contains two separate steps separated by $`U`$ instead of a single step broadened by $`U`$.

For a molecule chemically bonded to a metallic surface, e.g., a PDT molecule bonded by a thiol group to a gold surface, the broadening $`\mathrm{\Gamma }`$ is expected to be of the same magnitude or larger than $`U`$. This washes out CB effects as shown in Fig. 17. Therefore, the CB is not expected in this case. However, if the coupling to both contacts is weak we should keep the possibility of CB and the importance of unrestricted solutions in mind.

### 3.3 Broadening

So far we have treated the level $`ϵ`$ as discrete, ignoring the broadening $`\mathrm{\Gamma }=\mathrm{\Gamma }_1+\mathrm{\Gamma }_2`$ that accompanies the coupling to the contacts. To take this into account we need to replace the discrete level with a Lorentzian density of states $`D(E)`$:

$$D(E)=\frac{1}{2\pi }\frac{\mathrm{\Gamma }}{(Eϵ)^2+(\mathrm{\Gamma }/2)^2}$$
(16)
As we will see later, $`\mathrm{\Gamma }`$ is in general *energy-dependent* so that $`D(E)`$ can deviate significantly from a Lorentzian shape. We modify Eqs. 9, 10 for $`N`$ and $`I`$ to include an integration over energy:

$`N`$ $`=`$ $`2{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}\text{d}ED(E){\displaystyle \frac{\mathrm{\Gamma }_1f(E,\mu _1)+\mathrm{\Gamma }_2f(E,\mu _2)}{\mathrm{\Gamma }_1+\mathrm{\Gamma }_2}}`$ (17)
$`I`$ $`=`$ $`{\displaystyle \frac{2e}{\mathrm{}}}{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}\text{d}ED(E){\displaystyle \frac{\mathrm{\Gamma }_1\mathrm{\Gamma }_2}{\mathrm{\Gamma }_1+\mathrm{\Gamma }_2}}\left(f(E,\mu _1)f(E,\mu _2)\right)`$ (18)
The charging effect is included as before by letting the center $`ϵ`$, of the molecular density of states, float up or down according to Eqs. (12,13) for the restricted model or Eqs. (14,15) for the unrestricted model.

For the restricted model, the only effect of broadening is to smear out the I-V characteristics as evident from Fig. 15. The same is true for the unrestricted model as long as the broadening is much smaller than the charging energy (Fig. 16). But moderate amounts of broadening can destroy the Coulomb blockade effects completely and make the I-V characteristics look identical to the restricted model (Fig. 17). With this in mind, we will use the restricted model in the remainder of this chapter.

## 4 Non Equilibrium Green’s Function (NEGF) Formalism

The one-level toy model described in the last section includes the three basic factors that influence molecular conduction, namely, $`E_fϵ_0`$, $`\mathrm{\Gamma }_{1,2}`$ and $`U`$. However, real molecules typically have multiple levels that often broaden and overlap in energy. Note that the two-level model (Fig. 12) in the last section treated the two levels as *independent* and such models can be used only if the levels do not overlap. In general we need a formalism that can do justice to multiple levels with arbitrary broadening and overlap. The non-equilibrium Green’s function (NEGF) formalism described in this section does just that.

In the last section we obtained equations for the number of electrons, $`N`$ and the current, $`I`$ for a one-level model with broadening. It is useful to rewrite these equations in terms of the Green’s function $`G(E)`$ which is defined as follows:

$$G(E)=\left(Eϵ+i\frac{\mathrm{\Gamma }_1+\mathrm{\Gamma }_2}{2}\right)^1$$
(19)
The density of states $`D(E)`$ is proportional to the spectral function $`A(E)`$ defined as:

$`A(E)`$ $`=`$ $`2\text{Im}\left\{G(E)\right\}`$ (20)
$`D(E)`$ $`=`$ $`{\displaystyle \frac{A(E)}{2\pi }}`$ (21)
while the number of electrons, $`N`$ and the current, $`I`$ can be written as:

$`N`$ $`=`$ $`{\displaystyle \frac{2}{2\pi }}{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}\text{d}E\left(\left|G(E)\right|^2\mathrm{\Gamma }_1f(E,\mu _1)+\left|G(E)\right|^2\mathrm{\Gamma }_2f(E,\mu _2)\right)`$ (22)
$`I`$ $`=`$ $`{\displaystyle \frac{2e}{h}}{\displaystyle \underset{\mathrm{}}{\overset{\mathrm{}}{}}}\text{d}E\mathrm{\Gamma }_1\mathrm{\Gamma }_2\left|G(E)\right|^2\left(f(E,\mu _1)f(E,\mu _2)\right)`$ (23)
In the NEGF formalism the single energy level $`ϵ`$ is replaced by a Hamiltonian matrix $`\left[H\right]`$ while the broadening $`\mathrm{\Gamma }_{1,2}`$ is replaced by a complex energy-dependent self-energy matrix $`\left[\mathrm{\Sigma }_{1,2}(E)\right]`$ so that the Green’s function becomes a matrix given by:

$$G(E)=\left(ESH\mathrm{\Sigma }_1\mathrm{\Sigma }_2\right)^1$$
(24)
where $`S`$ is the identity matrix of the same size as the other matrices and the broadening matrices $`\mathrm{\Gamma }_{1,2}`$ are defined as the imaginary (more correctly as the anti-Hermitian) parts of $`\mathrm{\Sigma }_{1,2}`$:

$$\mathrm{\Gamma }_{1,2}=i\left(\mathrm{\Sigma }_{1,2}\mathrm{\Sigma }_{1,2}^{}\right)$$
(25)
The spectral function is the anti-Hermitian part of the Green’s function:

$$A(E)=i\left(G(E)G^{}(E)\right)$$
(26)
from which the density of states $`D(E)`$ can be calculated by taking the trace:

$$D(E)=\frac{\text{Tr}\left(AS\right)}{2\pi }$$
(27)
The density matrix $`\left[\rho \right]`$ is given by, c.f., Eq. 22:

$$\rho =\frac{1}{2\pi }\underset{\mathrm{}}{\overset{\mathrm{}}{}}\left[f(E,\mu _1)G\mathrm{\Gamma }_1G^{}+f(E,\mu _2)G\mathrm{\Gamma }_2G^{}\right]\text{d}E$$
(28)
from which the total number of electrons, N can be calculated by taking a trace:

$$N=\text{Tr}\left(\rho S\right)$$
(29)
The current is given by, c.f., Eq. 23:

$$I=\frac{2e}{h}\underset{\mathrm{}}{\overset{\mathrm{}}{}}\left[\text{Tr}\left(\mathrm{\Gamma }_1G\mathrm{\Gamma }_2G^{}\right)\left(f(E,\mu _1)f(E,\mu _2)\right)\right]\text{d}E$$
(30)
Equations 24 through 30 constitute the basic equations of the NEGF formalism which have to be solved self consistently with a suitable scheme to calculate the self-consistent potential matrix $`\left[U_{SCF}\right]`$, c.f., Eq. 13:

$$H=H_0+U_{SCF}$$
(31)
where $`H_0`$ is the bare Hamiltonian (like $`ϵ_0`$ in the toy model) and $`U_{SCF}`$ is an appropriate functional of the density matrix $`\rho `$:

$$U_{SCF}=F\left(\rho \right)$$
(32)
This self-consistent procedure is essentially the same as in Fig. 10 for the one level toy model, except that scalar quantities have been replaced by matrices:

$`ϵ_0`$ $``$ $`\left[H_0\right]`$ (33)
$`\mathrm{\Gamma }`$ $``$ $`\left[\mathrm{\Gamma }\right],\left[\mathrm{\Sigma }\right]`$ (34)
$`N`$ $``$ $`\left[\rho \right]`$ (35)
$`U_{SCF}`$ $``$ $`\left[U_{SCF}\right]`$ (36)
The sizes of all these matrices is ($`n\times n`$), $`n`$ being the number of basis functions used to describe the *molecule*. Even the self-energy matrices $`\mathrm{\Sigma }_{1,2}`$ are of this size although they represent the effect of infinitely large contacts. In the remainder of this section and the next section, we will describe the procedure used to evaluate the Hamiltonian matrix $`H`$, the self-energy matrices $`\mathrm{\Sigma }_{1,2}`$ and the functional “$`F`$” used to evaluate the self-consistent potential $`U_{SCF}`$ (see Eq. 32). But the point to note is that once we know how to evaluate these matrices, Eqs. 24 through 32 can be used straight forwardly to calculate the current.

*Non-orthogonal basis*: The matrices appearing above depend on the basis functions that we use. Many of the formulations in quantum chemistry use non-orthogonal basis functions and the matrix equations 24 through 32 are still valid as is, except that the elements of the matrix $`\left[S\right]`$ in Eq. 24 represents the overlap of the basis function $`\varphi _m(\overline{r})`$:

$$S_{mn}=\text{d}^3r\varphi _m^{}(\overline{r})\varphi _n(\overline{r})$$
(37)
For orthogonal bases, $`S_{mn}=\delta _{mn}`$ so that $`S`$ is the identity matrix as stated earlier. The fact that the matrix equations 24 through 32 are valid even in a non-orthogonal representation is not self-evident and is discussed in Ref. .

*Incoherent Scattering*: One last comment about the general formalism. The formalism as described above neglects all incoherent scattering processes inside the molecule. In this form it is essentially equivalent to the Landauer formalism . Indeed our expression for the current (Eq. 30) is exactly the same as in the transmission formalism with the transmission $`T`$ given by $`\text{Tr}\left(\mathrm{\Gamma }_1G\mathrm{\Gamma }_2G^{}\right)`$. But it should be noted that, the real power of the NEGF formalism lies in its ability to provide a first principles description of incoherent scattering processes - something we do not address in this chapter and leave for future work.

*A practical consideration*: Both Eq. 28 and 30 require an integral over all energy. This is not a problem in Eq. 30 because the integrand is non-zero only over a limited range where $`f(E,\mu _1)`$ differs significantly from $`f(E,\mu _2)`$. But in Eq. 28 the integrand is non-zero over a large energy range and often has sharp structures making it numerically challenging to evaluate the integral. One way to address this problem is to write:

$$\rho =\rho _{eq}+\mathrm{\Delta }\rho $$
(38)
where $`\rho _{eq}`$ is the equilibrium density matrix given by:

$$\rho =\frac{1}{2\pi }\underset{\mathrm{}}{\overset{\mathrm{}}{}}f(E,\mu )\left[G\mathrm{\Gamma }_1G^{}+G\mathrm{\Gamma }_2G^{}\right]\text{d}E$$
(39)
and $`\mathrm{\Delta }\rho `$ is the change in the density matrix under bias:

$$\rho =\frac{1}{2\pi }\underset{\mathrm{}}{\overset{\mathrm{}}{}}G\mathrm{\Gamma }_1G^{}\left[f(E,\mu _1)f(E,\mu )\right]+G\mathrm{\Gamma }_2G^{}\left[f(E,\mu _2)f(E,\mu )\right]\text{d}E$$
(40)
The integrand in Eq. 40 for $`\mathrm{\Delta }\rho `$ is non-zero only over a limited range (like Eq. 30 for $`I`$) and is evaluated relatively easily. The evaluation of $`\rho _{eq}`$ (Eq. 39) however still has the same problem but this integral (unlike the original Eq. 28) can be tackled by taking advantage of the method of contour integration as described in Ref. .

## 5 An Example: Quantum Point Contact (QPC)

Consider for example a gold wire stretched between two gold surfaces as shown in Fig. 18. One of the seminal results of mesoscopic physics is that such a wire has a quantized conductance equal to $`\frac{e^2}{\pi \mathrm{}}77.5\mu \text{A}/\text{V}\left(12.9\text{k}\mathrm{\Omega }\right)^1`$. This was first established using semiconductor structures  at $`4`$ K, but recent experiments on gold contacts have demonstrated it at room temperature . How can a wire have a resistance that is independent of its length? The answer is that this resistance is really associated with the interfaces between the narrow wire and the wide contacts. If there is scattering inside the wire it would give rise to an additional resistance in series with this fundamental interface resistance. The fact that a short wire has a resistance of $`12.9`$ k$`\mathrm{\Omega }`$ is a non-obvious result that was not known before 1988. This is a problem for which we do not really need a quantum transport formalism; a semi-classical treatment would suffice. The results we obtain here are not new or surprising. What is new is that we treat the gold wire as an $`Au_6`$ molecule and obtain well known results commonly obtained from a continuum treatment.

In order to apply the NEGF formalism from the last section to this problem, we need the Hamiltonian matrix $`[H]`$, the self energy matrices $`\mathrm{\Sigma }_{1,2}`$ and the self-consistent field $`U_{SCF}=F([\rho ])`$. Let us look at these one by one.

*Hamiltonian*: We will use a simple semi-empirical Hamiltonian which uses one s-orbital centered at each gold atom as the basis functions with the elements of the Hamiltonian matrix given by:

$`H_{ij}`$ $`=ϵ_0`$ $`\text{if}i=j`$
$`=t`$ $`\text{if}i,j\text{are nearest neighbors}`$
where $`ϵ_0=10.92`$ eV and $`t=2.653`$ eV. The orbitals are assumed to be orthogonal, so that the overlap matrix $`S`$ is the identity matrix.

*Self Energy*: Once we have a Hamiltonian for the entire molecule-contact system, the next step is to “partition” the device from the contacts and obtain the self-energy matrices $`\mathrm{\Sigma }_{1,2}`$ describing the effects of the contacts on the device. The contact will be assumed to be essentially unperturbed relative to the surface of a bulk metal so that the full Green’s function ($`G_T`$) can be written as (the energy $`E`$ is assumed to have an infinitesimal imaginary part $`i0^+`$):

$$G_T=\left(\begin{array}{cc}ESH& ES_{dc}H_{dc}\\ ES_{cd}H_{cd}& ES_cH_c\end{array}\right)^1=\left(\begin{array}{cc}G& G_{dc}\\ G_{cd}& G_C\end{array}\right)$$
(42)
where “$`c`$” denotes one of the contacts (the effect of the other contact can be obtained separately). We can use straight forward matrix algebra to show that:

$`G`$ $`=`$ $`\left(ESH\mathrm{\Sigma }\right)^1`$ (43)
$`\mathrm{\Sigma }`$ $`=`$ $`\left(ES_{dc}H_{dc}\right)\left(ES_cH_c\right)^1\left(ES_{cd}H_{cd}\right)`$ (44)
The matrices $`S_{dc}`$, $`S_c`$, $`H_{dc}`$, $`H_c`$ are all infinitely large since the contact is infinite. But the element of $`S_{dc}`$, $`H_{dc}`$ are non-zero only for a small number of contact atoms whose wavefunctions significantly overlap the device. Thus, we can write:

$$\mathrm{\Sigma }=\tau g_s\tau ^{}$$
(45)
where $`\tau `$ is the non-zero part of $`ES_{dc}H_{dc}`$ having dimensions ($`d\times s`$) where ’$`d`$’ is the size of the device matrix , and ’$`s`$’ is the number of surface atoms of the contact having a non-zero overlap with the device. $`g_s`$ is a matrix of size $`s\times s`$ which is a subset of the full infinite-sized contact Green’s function $`\left(ES_cH_c\right)^1`$. This surface Green’s function can be computed exactly by making use of the periodicity of the semi-infinite contact, using techniques that are standard in surface physics . For a one-dimensional lead, with a Hamiltonian given by Eq. 5, the result is easily derived :

$$g_s(E)=\frac{e^{ika}}{t}$$
(46)
where “$`ka`$” is related to the energy through the dispersion relation:

$$E=ϵ_02t\text{cos}(ka)$$
(47)
The results presented below were obtained using the more complicated surface Green’s function for an FCC (111) gold surface as described in Ref. . However, using the surface Green’s function in Eq. 46 gives almost identical results.

*Electrostatic potential*: Finally we need to identify the electrostatic potential across the device (Fig. 20) by solving the Poisson equation:

$$^2U_{tot}=\frac{e^2n}{ϵ_0}$$
(48)
with the boundary conditions given by the potential difference $`V_{app}`$ between the metallic contacts (here $`ϵ_0`$ is the dielectric constant). To simplify the calculations we divide the solution into a applied and self consistent potential ($`U_{tot}=U_{app}+U_{SCF}`$) where $`U_{app}`$ solves the Laplace equation with the known potential difference between the metallic contacts:

$`^2U_{app}=0`$ $`U_{app}=eV_n\text{on electrode ’n’}`$ (49)
Thus, $`U_{tot}`$ solves Eq. 48 if $`U_{SCF}`$ solves Eq. 48 with zero potential at the boundary.

$`^2U_{SCF}={\displaystyle \frac{e^2n}{ϵ_0}}`$ $`U_{SCF}=0\text{on all electrodes}`$ (50)
In the treatment of the electrostatic we assume the two contacts to be semi-infinite classical metals separated by a distance ($`W`$). This gives simple solutions to both $`U_{app}`$ and $`U_{SCF}`$. The applied potential is given by (capacitor):

$$U_{app}=\frac{V}{W}x$$
(51)
where x is the position relative to the midpoint between the contacts. The self consistent potential is easily calculated with the method of images where the potential is given by a sum over the point charges and all their images. However, to avoid the infinities associated with point charges, we adopt the Pariser-Parr-Pople (PPP) method  in the Hartree approximation. The PPP functional describing the electron-electron interactions is:

$$H_{ij}^{ee}=\delta _{ij}\underset{k}{}\left(\rho _{kk}\rho _{kk}^{eq}\right)\gamma _{ik}$$
(52)
where $`\rho `$ is the charge density matrix, $`\rho ^{eq}`$ the equilibrium charge density (in this case $`\rho _{ii}^{eq}=1`$ since we are modeling the s-electrons of gold) and the one center two-electron integral $`\gamma _{ij}`$. The diagonal elements $`\gamma _{ii}`$ are obtained from experimental data and the off-diagonal elements ($`\gamma _{ij}`$) are parameterized to describe a potential that decrease as the inverse of the distance ($`1/R_{ij}`$):

$$\gamma _{ij}=\frac{e^2}{4\pi ϵ_0R_{ij}+\frac{2e^2}{\gamma _{ii}+\gamma _{jj}}}$$
(53)
*Calculations on the QPC*: The results for the I-V and potential for a QPC are shown in Fig. 21. The geometry used was a linear chain of six gold atoms connected to the FCC (111) surface of the contacts in the center of a surface triangle. The Fermi energy of the isolated contacts was calculated to be $`E_f=8.67`$ eV, by requiring that there is one electron per unit cell.

As evident from the figure, the I-V characteristics is linear and the slope gives a conductance of $`77.3\mu \text{A}/\text{V}`$ close to the quantized value of $`\frac{e^2}{\pi \mathrm{}}77.5\mu \text{A}/\text{V}`$ as previously mentioned. What makes the QPC distinct from typical molecules is the strong coupling to the contacts which broadens all levels into a continuous density of states and any evidence of a conductance gap (Fig. 2) is completely lost. Examining the potential drop over the QPC shows a linear drop over the center of the QPC with slightly larger drop at the end atoms. This may seem surprising since transport is assumed to be “ballistic” and one expects no voltage drop across the chain of gold atoms. This can be shown to arise because the chain is very narrow (one atom in cross-section) compared to the screening length .

We can easily imagine a experimental situation where the device (QPC or molecule) is attached asymmetrically to the two contacts with one strong and one weak side. To model this situation we artificially decreased the interaction between the right contact and the QPC by a factor of $`0.2`$. The results of this calculation are shown in Fig. 22. A weaker coupling gives a smaller conductance, as compared with the previous figure. More interesting is that the potential drop over the QPC is asymmetric. Also in line with our classical intuition, the largest part of the voltage drop occurs at the weakly coupled contact with smaller drops over the QPC and at the strongly coupled contact. The consequences of asymmetric voltage drop over molecules has been discussed by Ghosh et al .

## 6 Concluding Remarks

In this chapter we have presented an intuitive description of the current-voltage (I-V) characteristics of molecules using simple toy models to illustrate the basic physics (sections 1-3). These toy models were also used to motivate the rigorous Non-Equilibrium Green’s Function (NEGF) theory (section 4). A simple example was then used in section 5 to illustrate the application of the NEGF formalism. The same basic approach can be used in conjunction with a more elaborate Hückel Hamiltonian or even an ab initio Hamiltonian. But for these advanced treatments we refer the reader to Refs. .

Some of these models are publicly available through the Purdue Simulation Hub (www.nanohub.purdue.edu) and can be run without any need for installation. In addition to the models discussed here, there is a Hückel model which is an improved version of the earlier model made available in 1999. Further improvements may be needed to take into account the role of inelastic scattering or polaronic effects, especially in longer molecules like DNA chains.

## Acknowledgments

Sections 2 and 3 are based on material from a forthcoming book by one of the authors (S.D.) entitled “Quantum Phenomena: From Atoms to Transistors”. It is a pleasure to acknowledge helpful feedback from Mark Ratner, Phil Bagwell and Mark Lundstrom. The authors are grateful to Prashant Damle and Avik Ghosh for helpful discussions regarding the ab initio models. This work was supported by the NSF under Grant No. 0085516 - EEC.

## Appendix A MATLAB Codes

The Matlab codes for the toy models can also be obtained at “www.nanohub.purdue.edu”.

### A.1 Discrete One Level Model

```
%  Toy model, one level
%  Inputs (all in eV)
E0=-5.5;Ef=-5;gam1=0.2;gam2=0.2;U=1;
%  Constants (all MKS, except energy which is in eV)
hbar=1.06e-34;q=1.6e-19;IE=(2*q*q)/hbar;kT=.025;
%  Bias (calculate 101 voltage points in [-4 4] range)
nV=101;VV=linspace(-4,4,nV);dV=VV(2)-VV(1);
N0=2/(1+exp((E0-Ef)/kT));
for iV=1:nV %  Voltage loop
  UU=0;dU=1;
  V=VV(iV);mu1=Ef-(V/2);mu2=Ef+(V/2);
  while dU>1e-6 %  SCF
    E=E0+UU;
    f1=1/(1+exp((E-mu1)/kT));f2=1/(1+exp((E-mu2)/kT));
    NN=2*((gam1*f1)+(gam2*f2))/(gam1+gam2); %  Charge
    Uold=UU;UU=Uold+(.05*((U*(NN-N0))-Uold));
    dU=abs(UU-Uold);[V UU dU];
  end
  curr=IE*gam1*gam2*(f2-f1)/(gam1+gam2);
  II(iV)=curr;N(iV)=NN;[V NN];
end
G=diff(II)/dV;GG=[G(1) G]; %  Conductance
h=plot(VV,II*10^6,’k’); %  Plot I-V

```

### A.2 Discrete Two Level Model

```
%  Toy model, two levels
%  Inputs (all in eV)
Ef=-5;E0=[-5.5 -1.5];gam1=[.2 .2];gam2=[.2 .2];U=1*[1 1;1 1];
% Constants (all MKS, except energy which is in eV)
hbar=1.06e-34;q=1.6e-19;IE=(2*q*q)/hbar;kT=.025;
n0=2./(1+exp((E0-Ef)./kT));
nV=101;VV=linspace(-6,6,nV);dV=VV(2)-VV(1);Usc=0;
for iV=1:nV
    dU=1;
    V=VV(iV);mu1=Ef+(V/2);mu2=Ef-(V/2);
    while dU>1e-6
        E=E0+Usc;
        f1=1./(1+exp((E-mu1)./kT));f2=1./(1+exp((E-mu2)./kT));
        n=2*(((gam1.*f1)+(gam2.*f2))./(gam1+gam2));
        curr=IE*gam1.*gam2.*(f1-f2)./(gam1+gam2);
        Uold=Usc;Usc=Uold+(.1*(((n-n0)*U’)-Uold));
        dU=abs(Usc-Uold);[V Usc dU];
    end
    II(iV)=sum(curr);N(iV,:)=n;
end
G=diff(II)/dV;GG=[G(1) G];
h=plot(VV,II); %  Plot I-V

```

### A.3 Broadened One Level Model

```
%  Toy model, restricted solution with broadening
%  Inputs (all in eV)
E0=-5.5;Ef=-5;gam1=0.2;gam2=0.2;U=1.0;
%  Constants (all MKS, except energy which is in eV)
hbar=1.06e-34;q=1.6e-19;IE=(2*q*q)/hbar;kT=.025;
%  Bias (calculate 101 voltage points in [-4 4] range)
nV=101;VV=linspace(-4,4,nV);dV=VV(2)-VV(1);
N0=2/(1+exp((E0-Ef)/kT));
for iV=1:nV %  Voltage loop
  UU=0;dU=1;
  V=VV(iV);mu1=Ef-(V/2);mu2=Ef+(V/2);
  nE=400; %  Nummerical integration over 200 points
  id=diag(eye(nE))’;
  EE=linspace(-10,0,nE);dE=EE(2)-EE(1);
  f1=1./(1+exp((EE-id*mu1)/kT));
  f2=1./(1+exp((EE-id*mu2)/kT));
  while dU>1e-4 %  SCF
    E=E0+UU;
    g=1./(EE-id*(E+i/2*(gam1+gam2)));
    NN=2*sum(g.*conj(g).*(gam1*f1+gam2*f2))/(2*pi)*dE;
    Uold=UU;UU=Uold+(.2*((U*(NN-N0))-Uold));
    dU=abs(UU-Uold);[V UU dU];
  end
  curr=IE*gam1*gam2*sum((f2-f1).*g.*conj(g))/(2*pi)*dE;
  II(iV)=real(curr);N(iV)=NN;[V NN curr E mu1 mu2];
end
G=diff(II)/dV;GG=[G(1) G]; %  Conductance
h=plot(VV,II,’.’); %  Plot I-V

```

### A.4 Unrestricted Discrete One-Level Model

```
% Toy model unrestricted solution
%  Inputs (all in eV)
E0=-5.5;Ef=-5;gam1=0.2;gam2=0.2;U=1;
%Constants (all MKS, except energy which is in eV)
hbar=1.06e-34;q=1.6e-19;IE=(q*q)/hbar;kT=.025;
%  Bias (calculate 101 voltage points in [-4 4] range)
nV=101;VV=linspace(-4,4,nV);dV=VV(2)-VV(1);
N0=1/(1+exp((E0-Ef)/kT));
for iV=1:nV %  Voltage loop
  U1=0;U2=1e-5;dU1=1;dU2=1;% Set U2=U1 for restricted solution
  V=VV(iV);mu1=Ef-(V/2);mu2=Ef+(V/2);
  while (dU1+dU2)>1e-6 %  SCF
    E1=E0+U1;E2=E0+U2;
    f11=1/(1+exp((E1-mu1)/kT));f21=1/(1+exp((E1-mu2)/kT));
    f12=1/(1+exp((E2-mu1)/kT));f22=1/(1+exp((E2-mu2)/kT));
    NN1=((gam1*f12)+(gam2*f22))/(gam1+gam2);
    NN2=((gam1*f11)+(gam2*f21))/(gam1+gam2);
    Uold1=U1;Uold2=U2;
    U1=Uold1+(.05*((2*U*(NN1-N0))-Uold1));
    U2=Uold2+(.05*((2*U*(NN2-N0))-Uold2));
    dU1=abs(U1-Uold1);dU2=abs(U2-Uold2);
  end
  curr1=IE*gam1*gam2*(f21-f11)/(gam1+gam2);
  curr2=IE*gam1*gam2*(f22-f12)/(gam1+gam2);
  I1(iV)=curr1;I2(iV)=curr2;
  N1(iV)=NN1;N2(iV)=NN2;[V NN1 NN2];
end
G=diff(I1+I2)/dV;GG=[G(1) G]; %  Conductance
h=plot(VV,I1+I2,’-’); %  Plot I-V

```

### A.5 Unrestricted Broadened One-Level Model

```
%  Toy model, unrestricted solution with broadening
%  Inputs (all in eV)
E0=-5.5;Ef=-5;gam1=0.2;gam2=0.2;U=1;
%  Constants (all MKS, except energy which is in eV)
hbar=1.06e-34;q=1.6e-19;IE=(q*q)/hbar;kT=.025;
%  Bias (calculate 101 voltage points in [-4 4] range)
nV=101;VV=linspace(-4,4,nV);dV=VV(2)-VV(1);
N0=1/(1+exp((E0-Ef)/kT));
nE=200; %  Nummerical integration over 200 points
id=diag(eye(nE))’;
EE=linspace(-9,-1,nE);dE=EE(2)-EE(1);
for iV=1:nV %  Voltage loop
  U1=0;U2=1;dU1=1;dU2=1;
  V=VV(iV);mu1=Ef-(V/2);mu2=Ef+(V/2);
  f1=1./(1+exp((EE-id*mu1)/kT));
  f2=1./(1+exp((EE-id*mu2)/kT));
  while (dU1+dU2)>1e-3 %  SCF
    E1=E0+U1;E2=E0+U2;
    g1=1./(EE-id*(E1+i/2*(gam1+gam2)));
    g2=1./(EE-id*(E2+i/2*(gam1+gam2)));
    NN1=sum(g1.*conj(g1).*(gam1*f1+gam2*f2))/(2*pi)*dE;
    NN2=sum(g2.*conj(g2).*(gam1*f1+gam2*f2))/(2*pi)*dE;
    Uold1=U1;Uold2=U2;
    U1=Uold1+(.2*((2*U*(NN2-N0))-Uold1));
    U2=Uold2+(.2*((2*U*(NN1-N0))-Uold2));
    dU1=abs(U1-2*U*(NN2-N0));dU2=abs(U2-2*U*(NN1-N0));
  end
  curr=IE*gam1*gam2*sum((f2-f1).*(g1.*conj(g1)+ ...
        g2.*conj(g2)))/(2*pi)*dE;
  II(iV)=real(curr);N(iV)=NN1+NN2;
  [V NN1 NN2 curr*1e6 E1 E2 mu1 mu2];
end
G=diff(II)/dV;GG=[G(1) G]; %  Conductance
h=plot(VV,II,’--’); %  Plot I-V

```


format: text

text:
# Counter-Intuitive Vacuum-Stimulated Raman Scattering

## Abstract

Vacuum-stimulated Raman scattering in strongly coupled atom-cavity systems allows one to generate free-running single photon pulses on demand. Most properties of the emitted photons are well defined, provided spontaneous emission processes do not contribute. Therefore, electronic excitation of the atom must not occur, which is assured for a system adiabatically following a dark state during the photon-generation process. We experimentally investigate the conditions that must be met for adiabatic following in a time-of-flight driven system, with atoms passing through a cavity and a pump beam oriented transverse to the cavity axis. From our results, we infer the optimal intensity and relative pump-beam position with respect to the cavity axis.

During the last years, many different types of single-photon sources have been proposed and successfully demonstrated Law96 ; Law97 ; Kim99 ; Lounis00 ; Kurtsiefer00 ; Brouri00 ; Michler00 ; Santori01 ; Yuan02 ; Michler00:2 ; Moreau01 ; Kuhn99 ; Hennrich00 ; Kuhn01 ; Kuhn02 . These activities were motivated by current quantum-cryptography schemes Tittel98 , which rely on the transmission of single photons between two parties to be secure, and by recent proposals Monroe02  on all-optical quantum information processing Knill01  and distributed quantum networking Cirac97 . For quantum networking, indistinguishable photons must be generated by a unitary, reversible process in order to transport single quantum bits from node to node. This last requirement is not met by most of the single-photon sources available so far. Only photon-generation schemes that are based on a deterministic and unitary energy exchange between a single emitter, e.g. an atom, and a single mode of the quantized radiation field Cirac97 ; Enk97:1 ; Parkins93 ; Parkins95 ; Pellizari95  might, in principle, meet the above requirements. Such a deterministic energy exchange has been successfully demonstrated in the microwave regime Rauschenbeutel99 ; Brattke01 , but these experiments do not provide free-running photons.

In this paper, we concentrate on a strongly coupled atom-cavity system in the optical domain, where the photons are transmitted through one of the cavity mirrors. The photons are generated by a vacuum-stimulated Raman transition Kuhn99 ; Hennrich00 ; Kuhn01 ; Kuhn02  that is adiabatically driven by a pump laser beam, interacting with the atoms while they are coupled to the cavity mode. The excitation scheme implements stimulated Raman scattering by adiabatic passage Vitanov01  (STIRAP), with the vacuum field of the cavity acting instead of a stimulating laser. Figure 1 illustrates this scheme: A cloud of Rubidium atoms in state $`|u`$, i.e. the $`5S_{1/2}(F=3)`$ hyperfine state of the electronic ground state, is released from a magneto-optical trap (MOT). The atoms fall through a 1 mm long cavity of Finesse $`F=\mathrm{60\hspace{0.17em}000}`$ which is placed 20 cm below the MOT, so that the atoms passing through the cavity have acquired a transverse (vertical) velocity of $`v=2`$m/s. One TEM00 mode of the cavity with waist $`w_C=35\mu `$m stimulates the Raman process, since it is resonant with the $`|e|g`$ transition of the atom, with $`|e`$ being an electronically excited state, and $`|g`$ being the other hyperfine ground state, $`5S_{1/2}(F=2)`$. The transition is pumped by a continuous laser beam of waist $`w_Pw_C`$, resonant with the $`|u|e`$ transition, that perpendicularly crosses the cavity axis and the trajectory of the falling atoms. With $`\delta _X`$ being the distance between the crossing points of cavity and beam axis, respectively, with the atomic trajectory, the Rabi frequency of the pump beam driving the $`|u|e`$ transition, $`\mathrm{\Omega }_P`$, and the vacuum-Rabi frequency induced by the cavity coupling the $`|e|g`$ transition, $`\mathrm{\Omega }_C`$, become time dependent and delayed due to the motion of the atoms, according to

$`\mathrm{\Omega }_P(t)`$ $`=`$ $`\mathrm{\Omega }_0e^{\left(\frac{vt+\delta _X}{w_P}\right)^2}\text{and}`$ (1)
$`\mathrm{\Omega }_C(t)`$ $`=`$ $`2g(t)=2g_0e^{\left(\frac{vt}{w_C}\right)^2},`$
where $`\mathrm{\Omega }_0`$ is the peak Rabi frequency of the pump beam and $`g_0`$ the peak atom-cavity coupling constant. In case of positive $`\delta _X`$, the atoms are first exposed to the pump beam, while the interaction with the cavity mode is delayed by $`\delta _X/v`$. The opposite situation is met in case of negative $`\delta _X`$. If we consider a Raman-resonant excitation with the detunings, $`\mathrm{\Delta }`$, of the cavity and the pump beam from their respective transitions being equal, the interaction Hamiltonian of the coupled laser-atom-cavity system reads (in a rotating frame)

$$H_{\mathrm{int}}=\frac{\mathrm{}}{2}\left[2\mathrm{\Delta }\sigma _{ee}2g(\sigma _{eg}a+a^{}\sigma _{ge})\mathrm{\Omega }_P(\sigma _{eu}+\sigma _{ue})\right].$$
(2)
Without coupling and for zero detuning, the three product states forming the one-photon manifold $`\{|u,0,|e,0,|g,1\}`$ are degenerate ($`|0`$ and $`|1`$ denote the photon number states of the cavity). In case of a coupled system, the degeneracy is lifted and the eigenfrequencies read

$$\omega ^0=0\text{and}\omega ^\pm =\frac{1}{2}\left(\mathrm{\Delta }\pm \sqrt{4g^2+\mathrm{\Omega }_P^2+\mathrm{\Delta }^2}\right).$$
(3)
The corresponding eigenstates are

$`|\varphi ^0`$ $`=`$ $`\mathrm{cos}\mathrm{\Theta }|u,0\mathrm{sin}\mathrm{\Theta }|g,1,`$ (4)
$`|\varphi ^+`$ $`=`$ $`\mathrm{cos}\mathrm{\Phi }\mathrm{sin}\mathrm{\Theta }|u,0\mathrm{sin}\mathrm{\Phi }|e,0+\mathrm{cos}\mathrm{\Phi }\mathrm{cos}\mathrm{\Theta }|g,1,`$
$`|\varphi ^{}`$ $`=`$ $`\mathrm{sin}\mathrm{\Phi }\mathrm{sin}\mathrm{\Theta }|u,0+\mathrm{cos}\mathrm{\Phi }|e,0+\mathrm{sin}\mathrm{\Phi }\mathrm{cos}\mathrm{\Theta }|g,1,`$
where the mixing angles $`\mathrm{\Theta }`$ and $`\mathrm{\Phi }`$ are given by

$$\mathrm{tan}\mathrm{\Theta }=\frac{\mathrm{\Omega }_P}{2g}\text{and}\mathrm{tan}\mathrm{\Phi }=\frac{\sqrt{4g^2+\mathrm{\Omega }_P^2}}{\sqrt{4g^2+\mathrm{\Omega }_P^2+\mathrm{\Delta }^2}\mathrm{\Delta }}.$$
(5)
Recently, we have proposed Kuhn99  and demonstrated Hennrich00 ; Kuhn01 ; Kuhn02  that the dark state $`|\varphi ^0`$ can be used to generate photons in the cavity in a deterministic way. To do so, the interaction of a single atom with the cavity must be strong when the interaction starts, so that

$$|u,0|\varphi ^0|=1\text{for}2g\mathrm{\Omega }_P.$$
(6)
Provided the system’s state vector, $`|\mathrm{\Psi }`$, adiabatically follows the dark state, $`|\varphi ^0`$, throughout the interaction, the contribution of $`|g,1`$ increases with rising pump Rabi frequency $`\mathrm{\Omega }_P`$, so that the transient photon-emission rate from the cavity reads

$$R_{emit}=2\kappa |g,1|\varphi ^0|^2\text{tr}(\rho )=2\kappa \frac{\mathrm{\Omega }_P^2}{4g^2+\mathrm{\Omega }_P^2}\text{tr}(\rho ),$$
(7)
where $`\kappa `$ is the field decay rate of the cavity and $`\rho =|\mathrm{\Psi }\mathrm{\Psi }|`$ the density matrix of the considered three-level system, which is not closed. Hence, $`\text{tr}(\rho )`$ is not one, but denotes the population remaining in the coupled system at any given moment. With the other eigenstates, $`|\varphi ^\pm `$, not being populated, transverse spontaneous emission losses from the excited state $`|e`$ do not occur, and the overall photon emission probability is

$$P_{emit}=_0^\tau 𝑑tR_{emit}1\text{for}2\kappa \tau \frac{\mathrm{\Omega }_P^2}{4g^2+\mathrm{\Omega }_P^2}1,$$
(8)
where $`\mathrm{}`$ denotes the average of its argument in the relevant time interval, $`[0\mathrm{}\tau ]`$. Obviously, the process can reach a photon-generation efficiency close to unity, provided adiabatic following is assured. Therefore, the set of parameters $`\{\mathrm{\Omega }_0,g_0,w_C,w_P,\mathrm{\Delta },\delta _X\}`$ must be chosen in such a way that the adiabaticity criterion Messiah59

$$|\omega ^0\omega ^\pm ||\varphi ^\pm |\frac{d}{dt}|\varphi ^0|=|\dot{\mathrm{\Theta }}|$$
(9)
is met throughout the whole process. This requires a careful adjustment of the displacement, $`\delta _X`$, between pump laser beam and cavity axis.

Figure 2 shows $`|(\omega ^0\omega ^\pm )/\dot{\mathrm{\Theta }}|`$ as a function of $`t`$ and $`\delta _X`$. The smaller this ratio, the higher the probability to have a non-adiabatic evolution. Obviously, for any choice of $`\delta _X`$, such situations are encountered. The evolution is always non-adiabatic long before and long after the atom interacts with laser beam and/or cavity. However, since the coupling constant and the Rabi frequency are both weak in these regimes and $`\mathrm{\Theta }`$ changes only slowly, the population in the atomic bare states is not affected. A very crucial situation is met when $`|\dot{\mathrm{\Theta }}|`$ reaches its maximum, i.e. when the atom crosses the dashed line in fig. 2, which lies between cavity axis and pump beam. For $`|\delta _X|>\mathrm{max}(w_C,w_P)`$, the level splitting $`|(\omega ^0\omega ^\pm )|`$ is too small to assure adiabaticity and losses to other eigenstates occur. Only in the intermediate regime, where the pump beam partially overlaps the cavity mode, i.e. for $`|\delta _X/2|\mathrm{min}(w_C,w_P)`$, adiabaticity is assured even for large $`|\dot{\mathrm{\Theta }}|`$, and the behaviour of the system is predetermined by the projection of the initial state $`|u,0`$ onto the eigenstates as soon as the adiabatic regime is entered. In this case, three major scenarios must be distinguished: (a) $`\delta _X>0`$. The atom first interacts with the pump beam, the initial state projects onto $`|\varphi ^\pm `$, and the photon is lost by transverse spontaneous emission. (b) $`\delta _X0`$. The dark state $`|\varphi ^0`$ is partially populated, so that the probability of a photon emission from the cavity might reach 50%. (c) $`\delta _X<0`$. The atom first interacts with the cavity, the initial state $`|u,0`$ projects onto the dark state $`|\varphi ^0`$, so that no losses occur except the desired photon emission from the cavity with a probability approaching 100%.

Figure 3 shows the number of photons emitted from the cavity as a function of pump beam displacement $`\delta _X`$ for five different pump intensities. The photons are counted during $`5.24`$ms while a single cloud of atoms passes through the cavity. Pump beam and cavity are Raman resonant, with a common detuning of $`\mathrm{\Delta }=10`$MHz from the direct atomic transitions. The pump beam displacement is adjusted using a piezoelectric mirror assembly, and beam position and waist are monitored by a CCD camera. The trace that belongs to the smallest pump intensity is used to calibrate the $`\delta _X`$-origin, assuming that the small Rabi frequency gives rise to an adiabatic evolution only for $`\delta _X=0`$, so that the peak photon number is found there. From this position, the peak emission shifts towards negative $`\delta _X`$ with increasing Rabi frequency, since the level splitting $`|\omega ^0\omega ^\pm |`$ and therefore the adiabatic regime both increase with $`\mathrm{\Omega }_0`$. This allows to pull pump beam and cavity further apart without loosing adiabaticity, so that the fraction of the initial state that projects onto the dark state $`|\varphi ^0`$ when the adiabatic regime is entered increases as well. The data also reveal that a significant overlap between pump beam and cavity mode is mandatory, otherwise adiabaticity is not assured when the atom reaches the point of maximum $`|\dot{\mathrm{\Theta }}|`$, which is located between cavity mode and pump beam axis, so that the number of emitted photons decreases again for $`\delta _X\mathrm{}`$. In our experiment, the optimum is found for $`\mathrm{\Omega }_Pg_0`$ and $`\delta _X(w_C+w_P)/2`$.

A comparison of the experimental results with a numerical simulation of the process must take into account that the atomic trajectories are not controlled. Therefore the effective atom-cavity coupling depends on the random ‘point of impact’, $`\stackrel{}{r}=(y,z)`$, of every single atom with respect to the cavity’s mode function, so that

$$g_0(y,z)g_0^{max}\mathrm{cos}(2\pi z/\lambda )\mathrm{exp}\left((y/w_C)^2\right).$$
(10)
The variation of $`w_C`$ along $`z`$ is neglected here, since the Rayleigh length of the cavity mode exceeds the cavity length. The dependence of the photon emission probability $`P_{emit}`$ on $`g_0`$ is highly non-linear, so that $`g_0`$ cannot be replaced by its mean value prior to a numerical calculation of $`P_{emit}`$. Instead,we have to average $`P_{emit}`$ over all possible points of impact, according to

$$\overline{P}_{emit}=_0^\lambda 𝑑z_{S_y/2}^{S_y/2}𝑑y\frac{P_{emit}(g_0(y,z))}{\lambda S_y}.$$
(11)
Note that $`S_y=100\mu `$m is the width of a slit aperture that is installed above the cavity. This aperture constrains the atomic trajectories perpendicular to the cavity axis to the interval $`[S_y/2\mathrm{}+S_y/2]`$. To calculate $`P_{emit}(g_0(y,z))`$, we numerically solve the master equation of the coupled system Kuhn99  and integrate (11). This calculation takes the cavity-field decay rate, $`\kappa `$, and the decay rate of the excited atomic level, $`\gamma `$, into account.

Figure 4 shows a simulation of the average $`\overline{P}_{emit}`$ as a function of $`\delta _X`$ for the range of pump intensities investigated in the experiment (see Fig. 3). For the smallest Rabi frequency the maximum photon emission probability is found at $`\delta _X=0`$, which justifies our $`\delta _X`$ calibration. Moreover, the simulation shows the same trend that is observed in the experiment, with the peak emission probability shifting towards negative $`\delta _X`$ with increasing Rabi frequency. However, simulation and experiment do not agree perfectly, and also the experimentally observed photon-number reduction at high pump intensities is not reproduced. These discrepancies lead to the conclusion that either the pump beam deviates from an ideal Gaussian beam in its wings, or that weak stray light of the pump beam (e.g. from the vacuum viewports) also hits the atoms. Both effects give rise to an electronic excitation of the atoms if the pump beam is very intense, which results in an early loss of photons by spontaneous emission into transverse modes.

From the theoretical considerations, the experimental results and the numerical calculations we draw the conclusion that vacuum-stimulated Raman transitions are most effectively driven by an adiabatic passage that results from a counter-intuitive interaction sequence, where the atoms are first coupled to the vacuum field of an empty cavity, stimulating the transition, and then exposed to a pump laser beam. To assure adiabaticity and to avoid losses to other states, a significant overlap of cavity mode and pump beam is required. Optimum conditions are found experimentally for a pump Rabi frequency that equals the maximum atom-cavity coupling constant, $`g_0^{max}`$, and for a beam displacement that equals the average waist of cavity and pump beam.

###### Acknowledgements.

This work was supported by the focused research program ‘Quantum Information Processing’ of the Deutsche Forschungsgemeinschaft, and by the European Union through the IST (QUBITS) and IHP (QUEST) programs.


text:
# Proof of the generalized Lieb-Wehrl conjecture for integer indices larger than one

## I Introduction

The Husimi function husimi  is a useful tool to investigate dynamical properties of quantum systems. Wehrl proposed a classical entropy wehrl  using the Husimi function and conjectured that the minimum of the entropy is obtained for coherent states. His conjecture was soon proved by Lieb lieb . However, Lieb’s analogous conjecture for spin coherent state in the same paper has remained unproved for more than 20 years.

The Husimi function for the density operator $`\widehat{\rho }`$ is defined as

$`\rho (\zeta )\zeta |\widehat{\rho }|\zeta `$ (1)
where $`|\zeta `$ is a coherent state. The generalized coherent states defined by Perelomov perelomov  can also be used instead of the standard coherent state. Gnutzmann and Życzkowski gnutzmann  have recently generalized the Wehrl entropy by analogy with the well-known Rényi entropy renyi . The Rényi-Wehrl entropy is defined as

$`S^{(q)}[\widehat{\rho }]={\displaystyle \frac{1}{1q}}\mathrm{log}M^{(q)}[\widehat{\rho }].`$ (2)
Here, $`M^{(q)}`$ is the moment of the Husimi function

$`M^{(q)}[\widehat{\rho }]=c{\displaystyle 𝑑\mu (\zeta )\{\rho (\zeta )\}^q}`$ (3)
where $`d\mu (\zeta )`$ is the Haar measure and $`c`$ is a normalization constant. In the limit $`q1`$, it reproduces the usual form of the entropy

$`\underset{q1}{lim}S^{(q)}[\widehat{\rho }]=S[\widehat{\rho }]c{\displaystyle 𝑑\mu (\zeta )\rho (\zeta )\mathrm{log}\rho (\zeta )}.`$ (4)
The generalized Lieb-Wehrl conjecture formed in gnutzmann  is the following.

###### Conjecture 1

The minimum value of the Rényi-Wehrl entropy for $`q>0`$ is obtained for coherent states.

We prove this conjecture for integer indices $`q2`$ in the cases of the coherent states of compact semisimple Lie groups. The case of spin ($`SU(2)`$) coherent states has already been proved gnutzmann ; schupp . In the case of the standard coherent state, $`q=1`$ corresponds to the original Wehrl’s conjecture proved in lieb , and theorem 3 therein gives the proof for arbitrary $`q>1`$.

## II Generalized coherent states

Before coming to the proof, let us review the definition and some properties of the generalized coherent states perelomov . Let $`G`$ be a compact semisimple Lie group. The Lie algebra of $`G`$ can be written in so-called Cartan basis $`\{H_j,E_\alpha \}`$. An irreducible representation of $`G`$ is characterized by the lowest (or the highest) weight. Let $`|\lambda `$ be the lowest weight state with the weight $`\lambda `$. We denote the irreducible representation space specified by $`|\lambda `$ as $`D_\lambda `$. Coherent states in $`D_\lambda `$ are obtained by the action of the group $`G`$ on $`|\lambda `$, which can be written explicitly as 111In fact, (5) cannot represent the coherent states orthogonal to $`|\lambda `$. However, it does not matter when we consider entropies such as (2) because the measure of the set of such states is zero.

$`|\zeta =𝒩(\zeta )\mathrm{exp}(\zeta _\alpha E_\alpha )|\lambda .`$ (5)
Here, $`𝒩(\zeta )`$ is the normalization constant and $`\alpha `$ runs over all positive roots. The expansion of the exponential function becomes a finite series because $`D_\lambda `$ is finite dimensional. Therefore

$`\stackrel{~}{\psi }(\zeta ){\displaystyle \frac{\psi |\zeta }{𝒩(\zeta )}}`$ (6)
is a polynomial of $`\zeta `$ for any state $`|\psi `$. We will use this fact later to prove lemma 1.

The “resolution of unity”

$`I_{D_\lambda }=\mathrm{dim}D_\lambda {\displaystyle 𝑑\mu (\zeta )|\zeta \zeta |}`$ (7)
is valid, and hence the normalization constant in (3) should be taken as

$`c_{D_\lambda }=\mathrm{dim}D_\lambda ,`$ (8)
in order to satisfy the normalization condition $`M^{(1)}=1`$.

## III Proof of the main result

Turning now to the proof, we first note the following lemma.

###### Lemma 1

If $`\widehat{\rho }`$ minimizes $`S^{(q)}[\widehat{\rho }]`$ for some $`q>0`$, $`\widehat{\rho }`$ must be a pure state.

This is a generalization of lemma 2 in lieb , and can be proved in the same way (see the appendix). Hereafter we concentrate on pure states.

Let us consider the simplest case $`q=2`$. We will see later that the discussion in this case can be easily generalized. We have to prove that the maximum of the second moment of the Husimi function

$`M_{|\phi }^{(2)}c_{D_\lambda }{\displaystyle 𝑑\mu (\zeta )|\zeta |\phi |^4}`$ (9)
is given by coherent states.

Let us rewrite (9) as

$`M_{|\phi }^{(2)}`$ $`=`$ $`\left(\phi |\phi |\right)c_{D_\lambda }{\displaystyle 𝑑\mu (\zeta )\left(|\zeta |\zeta \right)\left(\zeta |\zeta |\right)\left(|\phi |\phi \right)}.`$ (10)
The key observation is that the tensor product of the coherent states $`|\zeta |\zeta `$ is a coherent state in the irreducible representation space $`D_{2\lambda }`$ whose lowest weight state is $`|\lambda |\lambda `$. The same idea is used in schupp  for $`SU(2)`$. It can be shown explicitly as

$`(H_jI+IH_j)|\lambda |\lambda `$ $`=`$ $`2\lambda _j|\lambda |\lambda `$ (11)
$`(E_\alpha I+IE_\alpha )|\lambda |\lambda `$ $`=`$ $`0`$ (12)
$`|\zeta |\zeta `$ $`=`$ $`𝒩(\zeta )^2\mathrm{exp}[\zeta _\alpha (E_\alpha I+IE_\alpha )]|\lambda |\lambda `$ (13)
where $`\alpha `$ runs over positive roots. Therefore $`𝑑\mu (\zeta )\left(|\zeta |\zeta \right)\left(\zeta |\zeta |\right)`$ is proportional to the projection operator to $`D_{2\lambda }`$ in $`D_\lambda D_\lambda `$, i.e.

$`c_{D_\lambda }{\displaystyle 𝑑\mu (\zeta )\left(|\zeta |\zeta \right)\left(\zeta |\zeta |\right)}={\displaystyle \frac{\mathrm{dim}D_\lambda }{\mathrm{dim}D_{2\lambda }}}P_{D_{2\lambda }}.`$ (14)
Here the proportional constant is determined from (7). Hence

$`M_{|\phi }^{(2)}`$ $`=`$ $`{\displaystyle \frac{\mathrm{dim}D_\lambda }{\mathrm{dim}D_{2\lambda }}}\left(\phi |\phi |\right)P_{D_{2\lambda }}\left(|\phi |\phi \right)`$ (15)
and we have

$`|\phi |\phi D_{2\lambda }M_{|\phi }^{(2)}=M_{max}^{(2)}{\displaystyle \frac{\mathrm{dim}D_\lambda }{\mathrm{dim}D_{2\lambda }}}.`$ (16)
Then we should prove

$`|\phi :\text{coherent state}|\phi |\phi D_{2\lambda }.`$ (17)
LHS $``$ RHS is obvious, because $`|\phi |\phi `$ is obtained by the action of $`G`$ on $`|\lambda |\lambda `$ which belongs to $`D_{2\lambda }`$ by definition.

Then we prove LHS $``$ RHS. Let us consider the quadratic Casimir operator $`C_2`$ of the group $`G`$. By taking an orthonormal basis $`\{T_a\}`$ of the Lie algebra, it can be written as

$`C_2=𝑻^2{\displaystyle \underset{a}{}}T_a^2.`$ (18)
In the tensor product space $`D_\lambda D_\lambda `$, the quadratic Casimir operator is written as

$`C_2^{(2)}{\displaystyle \underset{a}{}}(T_a1+1T_a)^2=𝑻^21+1𝑻^2+2{\displaystyle \underset{a}{}}T_aT_a.`$ (19)
Therefore

$`(\phi |\phi |)C_2^{(2)}(|\phi |\phi )=2<𝑻^2>+2<𝑻>^2`$ (20)
where the symbol $`<>`$ denotes the expectation value for $`|\phi `$.

Here we note the lemma proved in delbourgo :

###### Lemma 2 (Delbourgo and Fox)

The minimum of the uncertainty $`(\mathrm{\Delta }T)^2<𝐓^2><𝐓>^2`$ is obtained for coherent states

$`(\mathrm{\Delta }T)^2\text{is minimum}|\phi :\text{coherent state}.`$ (21)
Since $`C_2=𝑻^2`$ is a constant in $`D_\lambda `$, $`<𝑻>^2`$ is maximum when $`|\phi `$ is a coherent state, and hence $`C_2^{(2)}`$ takes the maximum value at $`D_{2\lambda }`$ because it contains the tensor product of the coherent states. Therefore

$`|\phi |\phi D_{2\lambda }<𝑻>^2\text{ is maximum}|\phi :\text{coherent state}.`$ (22)
This completes the proof of (17).

The generalization to the higher moments is straightforward. By considering the tensor product space $`D_\lambda ^q\stackrel{q}{\stackrel{}{D_\lambda D_\lambda \mathrm{}D_\lambda }}`$, one can easily show

$`M_{|\phi }^{(q)}{\displaystyle \frac{\mathrm{dim}D_\lambda }{\mathrm{dim}D_{q\lambda }}}\left(\phi |^q\right)P_{D_{q\lambda }}\left(|\phi ^q\right)`$ (23)
and

$`\left(\phi |^q\right)C_2^{(q)}\left(|\phi ^q\right)=q<𝑻^2>+q(q1)<𝑻>^2`$ (24)
in obvious notation. From these two formulae we can conclude

$`M_{|\phi }^{(q)}=M_{max}^{(q)}|\phi ^qD_{q\lambda }|\phi :\text{coherent state}.`$ (25)
Combining this and Lemma 1, we obtain the main result of this paper:

###### Theorem 1

Let $`S^{(q)}`$ be the Rényi-Wehrl entropy with an integer index $`q2`$ of a compact semisimple Lie group. The minimum value of $`S^{(q)}`$ is obtained for coherent states

$`S^{(q)}[\widehat{\rho }]\text{is minimum}\widehat{\rho }=|\phi \phi |,|\phi :\text{coherent state}.`$ (26)
The minimum value $`S_{min}^{(q)}`$ is given by

$`S_{min}^{(q)}={\displaystyle \frac{1}{1q}}\mathrm{log}\left({\displaystyle \frac{\mathrm{dim}D_\lambda }{\mathrm{dim}D_{q\lambda }}}\right)`$ (27)
where $`D_\lambda `$ is the irreducible representation space in which $`S^{(q)}`$ is defined.

## IV Minimum values for some special cases

Finally, let us examine the general formula for the minima (27) in some special cases. Cases (i), (ii), and (iii) in the following can be used to investigate the dynamics of bosons, fermions and distinguishable particles, respectively 222$`U(N)`$, instead of $`SU(N)`$, is used in sugita2  for bosons and fermions. $`U(N)`$ is not semisimple. However, $`U(N)U(1)\times SU(N)`$, and the Abelian subgroup $`U(1)`$ corresponds to the total phase which is irrelevant for the Husimi functions. sugita2 .

* $`G=SU(N)`$, $`D_\lambda =[1^m]`$ (symmetric product of $`m`$ fundamental representations)

$`\mathrm{dim}D_{q\lambda }=\mathrm{dim}[1^{qm}]={\displaystyle \frac{(N+qm1)!}{(N1)!(qm)!}}`$ (28)
$`S_{min}^{(q)}={\displaystyle \frac{1}{1q}}\mathrm{log}{\displaystyle \frac{\mathrm{\Gamma }(N+m)\mathrm{\Gamma }(qm+1)}{\mathrm{\Gamma }(N+qm)\mathrm{\Gamma }(m+1)}}.`$ (29)
Assuming that (29) holds also for non-integer $`q`$, we obtain

$`\underset{q1}{lim}S_{min}^{(q)}=m\{\psi (N+m)\psi (m+1)\}`$ (30)
where $`\psi (z)\frac{d}{dz}\mathrm{log}\mathrm{\Gamma }(z)`$ is the digamma function. This result agrees with the value of $`S^{(1)}`$ for coherent states given in slomczynski , except for the shift by $`\mathrm{ln}N`$ due to another normalization of the coherent states.
* $`G=SU(N)`$, $`D_\lambda =[m]`$ (anti-symmetric product of $`m`$ fundamental representations)

$`dimD_{q\lambda }=\mathrm{dim}[m^q]`$ $`=`$ $`{\displaystyle \underset{j=0}{\overset{m1}{}}}{\displaystyle \frac{(N+qj1)!j!}{(Nj1)!(q+j)!}}`$ (31)
$`S_{min}^{(q)}={\displaystyle \frac{1}{1q}}{\displaystyle \underset{j=0}{\overset{m1}{}}}\mathrm{log}{\displaystyle \frac{\mathrm{\Gamma }(N+1j)\mathrm{\Gamma }(q+j+1)}{\mathrm{\Gamma }(N+qj)\mathrm{\Gamma }(j+2)}}.`$ (32)
In this case, the minimum value for $`q=1`$ is expected to be

$`\underset{q1}{lim}S_{min}^{(q)}={\displaystyle \underset{j=0}{\overset{m1}{}}}\{\psi (Nj+1)\psi (j+2)\}.`$ (33)
* $`G=\stackrel{m}{\stackrel{}{SU(N)\times \mathrm{}\times SU(N)}}`$, $`D_\lambda =[1]^m`$

$`\mathrm{dim}D_{q\lambda }=\mathrm{dim}[1^q]^m=\left({\displaystyle \frac{(N+q1)!}{(N1)!q!}}\right)^m`$ (34)
$`S_{min}^{(q)}={\displaystyle \frac{m}{1q}}\mathrm{log}{\displaystyle \frac{\mathrm{\Gamma }(N+1)\mathrm{\Gamma }(q+1)}{\mathrm{\Gamma }(N+q)}}.`$ (35)
The minimum value for $`q=1`$ is expected to be

$`\underset{q1}{lim}S_{min}^{(q)}=m\{\psi (N+1)\psi (2)\}.`$ (36)

###### Acknowledgements.

The author thanks Professor Gnutzmann for letting him know about reference schupp . He also thanks Professor Życzkowski for helpful discussions.

## Appendix A Proof of Lemma 1

Here we prove lemma 1: if $`\widehat{\rho }`$ minimizes $`S^{(q)}[\widehat{\rho }]`$ for some $`q>0`$, $`\widehat{\rho }`$ must be a pure state.

We consider three cases (i) $`0<q<1`$, (ii) $`q=1`$ and (iii) $`1<q`$ separately. $`\widehat{\rho }`$ minimizes $`M^{(q)}`$ and $`S`$ in cases (i) and (ii), respectively, and maximizes $`M^{(q)}`$ in case (iii). Let us decompose the density operator as $`\widehat{\rho }=_i\lambda _i\widehat{\rho }_i`$, where $`_i\lambda _i=1`$ and $`\widehat{\rho }_i|\psi _i\psi _i|`$ is a pure state. Since $`m^{(q)}(x)=x^q(0<q<1)`$ and $`s(x)=x\mathrm{log}x`$ are concave, and $`m^{(q)}(x)=x^q(q>1)`$ is convex,

$`M^{(q)}[\widehat{\rho }]`$ $``$ $`{\displaystyle \underset{i}{}}\lambda _iM^{(q)}[\widehat{\rho }_i](0<q<1)`$ (37)
$`S[\widehat{\rho }]`$ $``$ $`{\displaystyle \underset{i}{}}\lambda _iS[\widehat{\rho }_i]`$ (38)
$`M^{(q)}[\widehat{\rho }]`$ $``$ $`{\displaystyle \underset{i}{}}\lambda _iM^{(q)}[\widehat{\rho }_i](q>1).`$ (39)
The equality must hold in at least one of the three inequalities. Then $`\rho _i(\zeta )=\rho _j(\zeta )`$, and hence $`\left|\stackrel{~}{\psi _i}(\zeta )\right|=\left|\stackrel{~}{\psi _j}(\zeta )\right|`$ in the notation of (6), almost everywhere for all $`i,j`$. Since $`\stackrel{~}{\psi }_i`$ is a polynomial, $`\stackrel{~}{\psi }_i=a\stackrel{~}{\psi }_j`$, and hence $`|\psi _i=a^{}|\psi _j`$, with $`|a|=1`$. Therefore $`\widehat{\rho }_i=\widehat{\rho }_j`$ for all $`i,j`$, which means $`\widehat{\rho }`$ is a pure state.


text:
# Topological charge of Center Vortices

## 1 Introduction

Center vortices provide an appealing picture of confinement . When center vortices are removed from the Yang-Mills ensemble, the string tension  and the quark condensate  are lost. Therefore one should expect that center vortices can also explain spontaneous breaking of chiral symmetry. The mechanism of spontaneous breaking of chiral symmetry seems to be tied to the topological properties of the gauge fields. My talk is devoted to the topological charge of center vortices and is mainly based on .

In a mathematically idealized way center vortices can be defined in a gauge invariant way in $`D`$-dimensional space-time as $`(D2)`$-dimensional closed hypersurfaces (boundaries) $`\mathrm{\Sigma }`$ of electromagnetic flux which contribute a (non-trivial) center element $`Z`$ to the Wilson loop when they are non-trivially linked to the latter:

$`P\mathrm{exp}\left[{\displaystyle \underset{𝒞}{}}A\left(\mathrm{\Sigma }\right)\right]=Z^{L(𝒞,\mathrm{\Sigma })}.`$ (1)
Here $`L(𝒞,\mathrm{\Sigma })`$ denotes the linking-number between the loop $`𝒞`$ and the hypersurfaces $`\mathrm{\Sigma }`$. Below I will consider $`D=4`$.

## 2 Intersection of center vortex surfaces

Depending on the position of the vortex surface $`\mathrm{\Sigma }`$ in the 4-dimensional spacetime manifold the flux of the vortex can be electric or magnetic or both. Generically, the vortex surface $`\mathrm{\Sigma }`$ evolves in time and at a fixed time the center vortex represents a closed loop of magnetic flux, i.e. the $`B`$-field is tangential to the loop. There are also purely spatial vortex surfaces existing at a single time instant only. These vortex surface carry only electric flux being normal to the spatial surface $`\mathrm{\Sigma }`$. Obviously, a non-zero topological charge

$`\nu ={\displaystyle \frac{1}{16\pi ^2}}{\displaystyle \mathrm{d}^4x\mathrm{tr}(F\stackrel{~}{F})}={\displaystyle \frac{1}{4\pi ^2}}{\displaystyle \mathrm{d}^4x\stackrel{}{E}(x)\stackrel{}{B}(x)}`$
arises when a generic vortex patch (evolving in time) and a purely spatial vortex patch intersect111Note, by the Bianchi identity $`\stackrel{}{}\times \stackrel{}{E}=_t\stackrel{}{B}`$ a time-dependent $`\stackrel{}{B}`$-field generates an electric field $`\stackrel{}{E}`$, which is however, perpendicular to the $`\stackrel{}{B}`$-field. Hence a time-dependent magnetic flux loop alone will generically not generate a topological charge, unless non-parallel loop segments intersect (see below). Note, that the dual relation $`\stackrel{}{}\times \stackrel{}{B}=_t\stackrel{}{E}`$ (which is part of the Yang-Mills equation of motion) is usually not satisfied, so that a time-dependent $`\stackrel{}{E}`$-field not necessarily generates also a $`\stackrel{}{B}`$-field. It is therefore not surprising that the topological charge of a center vortex is given by its self-intersection number $`\nu =\frac{1}{4}I(\mathrm{\Sigma },\mathrm{\Sigma })`$.

The self-intersection number $`I(\mathrm{\Sigma },\mathrm{\Sigma })`$ receives contributions from two types of singular points: (i) Transversal intersection points arising from the intersection of two different surface patches (see fig. 1) and (ii) twisting points occurring on a single surface patch twisting around a point in such a way to produce 4-linearly independent tangent vectors (see fig. 1). Transversal intersection points yield a contribution $`\pm 2`$ to the oriented intersection number, where the sign depends on the relative orientation of the two intersecting surface pieces. Twisting points yield always contributions of module smaller than 2.

Fig. 2 shows a vortex with a transversal intersection point $`\left(\nu =\frac{1}{2}\right)`$ and two twisting points $`\left(\nu =\frac{1}{8}\right)`$ at the front and back, respectively, edges of the configuration at the intermediate time $`(n_0=2)`$. Further twisting points $`\left(\nu =\frac{1}{8}\right)`$ occur at the initial $`\left(n_0=1\right)`$ and final $`\left(n_0=3\right)`$ times, so that the total topological charge of this configuration vanishes.

## 3 Writhing of center vortex loops

For generic center vortices representing closed magnetic flux loops $`𝒞(t)`$ evolving in time the topological charge can be expressed as

$`\nu ={\displaystyle \frac{1}{4}}{\displaystyle dt_tW\left(𝒞(t)\right)},`$ (2)
where $`W(𝒞)=L(𝒞,𝒞)`$ is the writhing number defined by the Gaussian linking number $`L(𝒞_1,𝒞_2)`$ between two loops $`𝒞_1`$ and $`𝒞_2`$.

The writhing number $`W\left(𝒞\right)`$ is a continuous function of the shape of the loop $`𝒞`$. It is not a topological invariant and thus not integer valued. It vanishes for planar curves or for curves possessing a symmetry plane, and in this sense measures the chirality of the loop. Furthermore, it suffers a discontinuity when two non-parallel line-segments of the curve cross.

Fig. 3 shows the time evolution of a closed magnetic vortex loop in ordinary 3-space, which on a 4-dimensional lattice gives rise to the configuration shown in fig. 2, after eliminating purely spatial vortex patches, which can be considered as lattice artifacts due to the discretization of time. For simplicity I have kept the cubistic form of the vortex in $`D=3`$ space, so that the loops consist of straight line segments.

Singular changes of the vortex loop occur at the creation and annihilation of the loop (which correspond to the twisting points at $`n_0=1,3`$ in fig. 2), where the writhing number changes by $`\mathrm{\Delta }W=\frac{1}{2}`$, and at the intermediate time $`t=\overline{t}_2`$, where the two long line-segments cross. This crossing changes $`W`$ by $`(2)`$ (i.e. $`\mathrm{\Delta }\nu =\frac{1}{2}`$) and corresponds in $`D=4`$ to the transversal intersection point at $`n_0=2`$, see fig. 2. Furthermore, when the two long loop segments cross at $`t=\overline{t}_2`$ the two short horizontal loop segments at the front and back edges reverse their direction, which can be interpreted as twisting these loop segments by an angle $`\pi `$. In the $`D=4`$ lattice realization of this vortex shown in fig. 2 these twistings correspond to the twisting points at $`n_0=2`$ at the front and back edges of the vortex.

It turns out that transversal intersection points do not give rise to a change of the so-called twist  of the vortex loop, while twisting points do, which justifies their name .

In view of the fact that center vortices carry spots of topological charge restricted to $`|\nu |\frac{1}{2}`$ the recent lattice measurement of the topological charge distribution  gives further support for the vortex picture of the QCD vacuum.

## 4 Quark guides

A non-vanishing total topological charge requires non-oriented vortex surfaces, which carry magnetic monopole loops at the boundary between oppositely oriented vortex patches . By the Athiya-Singer index theorem, $`\nu =N_LN_R`$, a non-zero topological charge $`\nu `$ is connected to the difference between the numbers $`N_{L/R}`$ of left and right handed quark zero-modes. Fig. 4 shows the probability density of the zero-modes of the quarks moving in the background of two pairs of intersecting center vortices on the 4-dimensional torus . As one observes the quark zero-modes are concentrated on the center vortex sheets and are in particular localized at the intersection points, the spots of topological charge $`\nu =\frac{1}{2}`$. If the quark zero-modes dominate the quark propagator, at low energies the quarks will travel along the center vortex sheets and can move from one vortex to an other through the intersection points. Since the center vortices percolate in the QCD vacuum, we expect also the percolation of the quark trajectories, which will eventually result in a condensation of the quarks.


text:
# A New Paradigm in Chaotic Ray Scattering

## I Introduction

Ray-chaos (exponential separation of nearby-originating rays, due to the geometrical features of the reflecting boundaries, and/or the refractive properties of the medium) is a lively subject, as witnessed by an increasing number of topical papers.

Well known paradigms of two-dimensional (henceforth 2D) ray-chaos include the Sinai Sinai\_bill  and Bunimovic Buni\_bill  billiards (internal boundary value problems) and the $`n`$-disk ($`n3`$) pinball-scatterer (external boundary value problem) pinball . These systems have been extensively studied both theoretically Theory\_Review , in view of their relevance in the context of quantum chaology Quantum\_chaology , and experimentally (see e.g. Stockmann , and Smilanski ), thanks to the formal equivalence between the 2D Schrödinger equation for a point particle in an infinite-wall potential, and the 2D (scalar) Helmholtz equation with perfectly conducting boundaries of Electromagnetics, which makes an experimental approach based on microwave measurements viable.

In this paper we introduce and analyze a new intriguing paradigm in ray-chaotic electromagnetic (henceforth EM) boundary-value problems. In its simplest (2D) form, the problem’s geometry is sketched in fig. 1, and consists of a dielectric layer with exponentially tapered refraction index matched at $`z=0`$ to the outer vacuum (viz., $`n(0)=1`$), and backed by a smooth perfectly-conducting periodic (undulated) surface $`z=\zeta (x)`$.

This geometry bears several elements of novelty. It is, perhaps, the simplest (2D) internal/external electromagnetic ray-chaotic boundary value problem featuring both guided and scattered rays at the same time. It also includes a (linear, isotropic, time-invariant, non-dispersive) dielectric layer as a basic ingredient at variance of billiards and pinballs, which are pure perfect-conductor-and-vacuum problems.

The 2D ray trajectories bear some similarity to the space-time ones of the well known (one-dimensional) Fermi-Pustil’nikov billiard Fermi\_bill , In this analogy, the gradient of the refractive-index plays the role of gravity, bending the ray-paths downward in the dielectric layer, while the undulated surface plays the role of the vibrating table, changing the slope of the ray at each bounce.

## II Ray Tracing

The system in Fig. 1 can be parametrized in terms of the dimensionless variables $`n_f=n(h)`$, $`a/h`$, $`\mathrm{\Delta }/a`$. In this section we illustrate the main relevant results of straightforward ray-tracing simulations run for the geometry of Fig. 1, with $`a/h=1`$ and $`\mathrm{\Delta }/a=2\pi /5`$.

In Fig. 2 three examples are shown of (multi-hop) ray-paths originating from rays with close-by incidence points ($`\mathrm{\Delta }x10^3a`$) and the same incidence angle ($`\theta _i=\pi /12`$). In all three cases, a rapidly increasing ray separation is observed, resulting into completely different exit-angles. It will be shown soon that the ray separation increases exponentially in time, as the rays bounce back and forth through the dielectric layer.

In Fig. 3 the exit-angle is shown vs. the (scaled) incidence point position $`\overline{x}=x/a`$, $`a`$ being the spatial period of the undulated surface, for a fixed incidence angle $`\theta _i=\pi /12`$, and a few representative values of $`n_f`$. Intervals of distinctly irregular behaviour are observed. The measure of the set of incidence points corresponding to regular behaviour is seen to depend on the value of $`n_f`$, and shrinks to zero as $`n_f`$ is increased. On the other hand, at fixed $`n_f`$, the measure of the set of incidence points corresponding to regular behaviour is found to be essentially independent of the incidence angle, as exemplified in Fig. 4, where $`n_f=15`$ and $`\theta _i=k\pi /12`$, $`k=1,3,5`$.

Closer investigation reveals that smaller and smaller intervals of regular behaviour are hidden within arbitrarily small intervals of $`\overline{x}`$, as shown in Fig. 5 for $`\theta _i=\pi /12`$.

As anticipated, in the irregular regime, the separation between neighbouring incident rays increases exponentially, as the rays bounce back and forth through the dielectric layer. This is shown in Fig. 6, where the logarithm of the ray-separation (in units of its initial value) is displayed as a function of time (sec) for $`n_f=15`$, and $`\theta _i=k\pi /12`$, $`k=1,3,5`$.

## III Full Wave Analysis: Results

A full-wave (Floquet mode) analysis (which is omitted here for brevity, and will be found in a forthcoming paper) suggests that in the limit of increasingly large $`a/\lambda _0`$, $`\lambda _0`$ being the field wavelength in vacuum both the internal $`(z0)`$ and the external ($`z>0`$) fields eventually comply with Berry’s conjecture Berry\_conj .

This is illustrated in Fig.s 7 and 8, which refer to the special (scalar) case where the electric field has only a nonzero component orthogonal to the $`xz`$ plane, and all other parameters hold the same values as in Sect. 2.

The field scattered into the half-space $`z>0`$ is indeed a superposition of (propagating or evanescent) plane waves whose wave-vectors span the characteristic (grating-lobe) directions Beck\_Spizz

$$\mathrm{sin}\theta _p=\mathrm{sin}\theta _in_f\frac{p\lambda _0}{a},p𝒵.$$
(1)
with (complex) amplitudes $`\mathrm{\Gamma }_pE_i`$ where $`E_i`$ is the (complex) amplitude of the incident plane wave, and $`\theta _i`$ is the incidence angle. It is seen from (1) that the number of propagating scattered plane waves is $`n_fa/\lambda _0`$. Numerical simulations further show that upon increasing $`a/\lambda _0`$, the distribution of the phase among the pertinent (complex) coefficients $`\mathrm{\Gamma }_p`$ becomes more and more evenly distributed in $`(0,2\pi )`$.

On the other hand, as seen from Fig. 7, the agreement between the (averaged average , normalized) spatial autocorrelation function $`E[\stackrel{}{r}+\widehat{x}(s/2)]E^{}[\stackrel{}{r}\widehat{x}(s/2)]`$ of the internal ($`z0`$) field and the $`J_0(n(z)k_0|\stackrel{}{s}|)`$ Bessel function limit-form predicted by Berry’s conjecture is remarkably good already at $`a/\lambda _0=6`$.

It is also noted that the complex coefficients $`\mathrm{\Gamma }_p`$ show a strong sensitivity to the incidence angle $`\theta _i`$, as illustrated in Fig. 8 for $`a/\lambda _0=9`$. In Fig. 8 the case of a plain sinusoidal surface is also shown ($`n_f=1`$), for which the $`\mathrm{\Gamma }_p`$s are seen to depend fairly smoothly on the incidence angle. This peculiar behaviour has a special notable consequence, also suggested by numerical simulations. In the limit where $`a/\lambda _00`$ the spread of the distribution of field values at any given point corresponding to a large number of different incidence angles in a sector of aperture $`\mathrm{\Delta }\theta _i`$. also remains $`𝒪(\mathrm{\Delta }\theta _i)`$, and the mean changes with the chosen field point. On the other hand, in the opposite limit where $`a/\lambda _0\mathrm{}`$, for any nonzero $`\mathrm{\Delta }\theta _i`$, albeit small, the spread of the distribution blows-up to a finite value, while both the spread and the mean of the distribution become almost independent from the chosen field point. This scenario is likely to occur in all wave-scattering problems exhibiting chaos in the $`\lambda 0`$ (ray) limit.

## IV Conclusions - Application Potential

Far from being a mere exercise of pure academical interest, the proposed scattering system holds the potential for several applications of practical interest.

The findings presented above suggest that an incident narrow coherent beam of electromagnetic radiation would be scattered off the system incoherently and isotropically if $`a/\lambda _01`$. Thus the proposed system might be used as a coating for radar-puzzling applications.

The multi-hop nature of the ray-paths also suggest that upon addition of suitably small lossy  dielectric and/or backing-surface losses the system might act as a smart radio-frequency absorber, featuring both fairly uniform power deposition in the dielectric bulk, and nicely small power reflection-coefficient for a wide interval of incidence angles.

## Acknowledgements

We acknowledge several stimulating discussions with professor Leo B. Felsen (Boston University, Boston MA, USA) and professor Giorgio Franceschetti (University of Naples, IT, and UCLA, Los Angeles CA, USA). This paper is dedicated to the memory of dr. Maurizio Mangrella.

Captions to the Figures

Fig 1 - Sketch of problem’s geometry. A perfectly conducting undulated surface with peak-to-peak height $`\mathrm{\Delta }`$ and spatial period $`a`$ is topped by a dielectric layer with maximum thickness $`h`$. The refractive index of the layer matches that of the vacuum at $`z=0`$ and grows exponentially with depth up to a value $`n_f`$ at $`z=h`$.

Fig. 2 - Multi-hop ray-paths originating from close-by incidence points display rapidly increasing separations, and eventually emerge with widely different exit angles.

Fig. 3 - Exit angle (deg.) vs. scaled incidence point position $`\overline{x}=x/a`$ for fixed incidence angle ($`\theta _i=\pi /12`$) and different values of $`n_f`$: $`n_f=1.4\text{(top)}`$, $`n_f=30\text{(mid)}`$, $`n_f=90\text{(bottom)}`$.

Fig. 4 - Exit angle (deg.) vs. scaled incidence point position $`\overline{x}=x/a`$ for fixed $`n_f`$ and different incidence angles. $`\theta _i=\pi /12\text{(top)}`$, $`\theta _i=3\pi /12\text{(mid)}`$, $`\theta _i=5\pi /12\text{(bottom)}`$.

Fig. 5 - Exit angle (deg.) vs. scaled incidence point position $`\overline{x}=x/a`$ for $`n_f=15`$ and $`\theta _i=\pi /12`$. Intervals corresponding to regular behaviour are observed at smaller and smaller scales (top to bottom).

Fig. 6 - Natural logarithm of separation between close-by incident rays as a function of time for $`n_f=15`$ and different incidence angles: $`\theta _i=\pi /12\text{(left)}`$, $`\theta _i=3\pi /12\text{(mid)}`$, $`\theta _i=5\pi /12\text{(right)}`$.

Fig. 7 - Full wave analysis with $`a/\lambda _0=6`$. Averaged normalized spatial autocorrelation function of the internal $`(z0)`$ field.

Fig. 8 - Full wave analysis with $`a/\lambda _0=9`$. Top: modulus (left) and phase (right) of $`\mathrm{\Gamma }_p`$ for two close-by incidence angles. Bottom: the same for a plain ($`n_f=1`$) undulated surface.


text:
# Quantum Computation with Vibrationally Excited Molecules

## Figure captions