Molecular reorganization energy in quantum-dot cellular automata switching

We examine the impact of the intrinsic molecular reorganization energy on switching in two-state quantum-dot cellular automata (QCA) cells. Switching a bit involves an electron transferring between charge centers within the molecule. This in turn causes the other atoms in the molecule to rearrange their positions in response. We capture this in a model that treats the electron motion quantum-mechanically, but the motion of nuclei semiclassically. This results in a non-linear Hamiltonian for the electron system. Interaction with a thermal environment is included by solving the Lindblad equation for the time-dependent density matrix. The calculated response of a molecule to the local electric field shows hysteresis during switching when the sweep direction is reversed. The relaxation of neighboring nuclei increases localization of the electron, which provides an intrinsic source of enhanced bistability and single-molecule memory. This comes at the cost of increased power dissipation.


I. INTRODUCTION
Quantum-dot cellular automata (QCA) is a potential alternative to traditional complimentary metal-oxide semiconductor (CMOS) transistor logic for nanoelectronics.2][3][4][5][6][7][8][9] A single mobile electron can be localized in a quantum dot, which is simply region of space with a surrounding potential barrier that quantizes the charge enclosed.The arrangement of charge among multiple dots in a QCA cell encodes information and electron transfer (ET) between dots through quantum mechanical tunneling enables switching.The electric field from one cell influences neighboring cells enabling device operation.A two-dot cell captures the most basic elements of QCA operation.Three-dot or six-dot cells can support clocked control of the flow of information, and allow power gain by replacing energy lost to dissipative processes through work done by the clock.
The QCA approach has been implemented in small metallic dots. 10Logic gates including majority gates, 11 digital latches, 12 and shift registers 13,14 have been demonstrated.Power gain 15 and signal restoration have also been demonstrated in these systems.
Semiconductor QCA cells have been demonstrated in GaAs 16 and silicon, 17 using dots formed by surface depletion gates.Another means of forming QCA cells in silicon is using a few implanted donor atoms to form the dots. 18olecular-scale QCA operation at room temperature has been achieved using dots formed by STM patterning of single dangling bonds on a hydrogen-passivated silicon surface. 19chieving QCA operation using single molecules [20][21][22] has focused attention on the richly-studied field of mixed-valence chemistry.In these molecules, the role of dots is played by charge centers, frequently formed by metal atoms surrounded by coordinating ligands.The central metal atom can exist in more than one charge state and can be reversibly oxidized a) Electronic mail: lent@nd.eduor reduced without breaking any chemical bonds.A second charge center within the molecule is connected through a bridging ligand so an electron can tunnel between dots.Of interest for QCA operation is the case when the charge centers are identical.If the bridging ligand presents a sufficient barrier to tunneling, the mobile charge will localize on one of the dots.The localization can be enhanced by the relaxation of the other atoms in molecule in response to the presence or absence of the charge on one or the other dot.
A simple example of such a mixed-valence molecule is the diferrocenylacetylene (DFA) 23 molecule shown in Figure 1.In DFA, a pair of ferrocenyl groups are connected through a bridging ligand.The two iron centers act as the two dots where an electron can be localized.
The potential advantages of single-molecule devices of this type are several.The functional density could be enormous because the footprint of each molecule on the surface is of the order of a square nanometer.In reality, it would not be practical to separately address and control each molecule separately, so clocking schemes and appropriate computational architectures, such as a Kogge network, 24 would be necessary.Another advantage is that the intrinsic speed of such devices could be very high.Electron transfer in mixed-valence molecules can occur at picosecond timescales. 25Finally, very low power dissipation may be possible because switching the FIG. 2. Schematic representation of effect of reorganization energy λ .A molecular double-dot can be viewed as a double potential well holding one mobile electron.If the electron is on the left (right) dot the polarization is +1 (-1).Electron occupation of a dot induces a relaxation of the positions of the surrounding atoms in response.This relaxation lowers the energy of the occupied dot compared to the unoccupied dot by the reorganization energy λ .bit requires only the motion of a single electron within the molecule.In contrast to molecular diodes or transistors, no current flow through the molecule is necessary.
Several candidate molecules for QCA have been synthesized and characterized.Fehlner et al. have synthesized QCA double-dots which were then covalently bonded to a semiconductor substrate. 26,27They demonstrated controlled switching of the molecular charge between dots using external electrodes.][30][31] Most mixed-valence species are ions.This may present complications due to the presence of neutralizing counterions that could bias the charge configuration of the QCA molecule.One strategy is to include an electron donor (or acceptor) at a symmetric position within the molecule.3][34] Christie et al. have synthesized such a molecule, which could form the basis of a clockable three-dot QCA cell. 35olecular QCA operation depends fundamentally on electric field driven ET between two dots within the molecule.We focus here on a simple model system that captures the essential features of this operation.Electron transfer in these systems is driven by an electric field, either produced by neighboring QCA molecules, or by clocking electrodes, or by input electrodes at the edges of a circuit array.This is in contrast to the usual situation in mixed valence chemistry, where ET is produced by either the random motion of solvent molecules in solution or by incident light.
Full ab initio quantum chemistry calculations of molecular states are possible for the static problem in an electric field, 7,[36][37][38][39] but calculation of time-dependent switching dynamics are too costly and one must use reduced-state models that capture the relevant physics.The parameters of the simplified system can be taken from either quantum chemistry calculations or from experiments.
The key parameters for such a mixed-valence double dot system model are: the inter-dot Hamiltonian tunneling matrix element, here called γ and often called H ab in the chemical literature; the on-site energy difference due to the applied field, here called ∆; the Marcus reorganization energy λ that captures the effect of the nuclear relaxation due to charge oc-cupancy of one or the other dot; the energy relaxation time T d which characterizes the energetic coupling between the molecule and the surrounding environment; and the environmental thermal energy k B T .
The inter-dot tunneling energy γ determines the intrinsic time scale for the system and can be varied over 12 orders of magnitude by adjusting the bridging ligand. 25,40Reference 31 shows the dramatic effect on localization of simply altering the geometry of the connection between the dots and the bridging cyclopentadiene ring from the meta to para configuration.
A molecular double-dot can be simply modelled as a double potential well along the spatial axis connecting the two metal centers.(Note that this is not the Marcus double well as a function of the reaction coordinate).The reorganization energy λ lowers the energy of the occupied well compared to the unoccupied well due to the relaxation of other atoms, either those in the molecule itself or in the surrounding environment, as shown schematically in Figure 2. The value of λ in solution can be affected by changing the solvent or by altering the chemical structure of the dot ligands.The reorganization energy can also be strongly affected by the attachment of the molecule to a surface.If the effective stiffness of the dot ligands is increased because of the constraint of the surface, the reorganization energy becomes smaller.(For a given force, a spring with a larger spring constant will store less energy than one with a smaller spring constant).A recent examination of the reorganization energy of a molecule with two ferrocene redox centers connected by a naphthalene linker and deposited on a NaCl substrate showed a reorganization energy in the few meV range, rather than the few hundred meV range for ferrocene in solution. 41,42In a similar way, the strength of the coupling to the thermal environment which determines T d is dependent on the details of the surface attachment.
The reorganization energy λ can be obtained theoretically by high quality quantum chemistry calculations.The nuclear motion of all the atoms in the molecule that are coupled to the electron transfer can be expressed in terms of the normal modes of vibration.These can then be treated quantum mechanically, though again this is computationally costly because of the need to account for the harmonic oscillator degrees of freedom. 43,44e employ a parameterized Hamiltonian to describe the electron moving between the two charge centers under a driving field.The electron motion is linearly coupled to the nuclear reorganization via a single vibrational mode (the antisymmetric breathing mode).The energy transferred to the nuclear degrees of freedom is accounted for semiclassically so energy is conserved in the electron-vibron system.This leads to a nonlinear Hamiltonian for the electron degrees of freedom.Finally, the interaction with the thermal environment is included through Lindblad operators.Our focus is on accounting for the "self-trapping" aspect of the relaxation of nuclear coordinates associated with ET, the associated hysteresis in switching behavior, and the effect on power dissipation.
Low power dissipation for device switching is a key requirement for any nanoelectronics.For classical systems, the dissipated energy exhibits an inverse linear dependence on FIG. 3. Schematic of a two-state quantum-dot system.The black circles represent the quantum dots, the lines connecting the dots represent the inter-dot tunneling paths, and the red circles denote the localized charge.When charge is localized on the left dot (L), the system is in quantum state |L and the polarization P = 1 representing a binary 0. When the charge is localized on the right dot (R), the system is in quantum state |R and the polarization P = −1 representing a binary 1.
switching time (1/T s ).To calculate the energy dissipated for a quantum system, a first-order approach is to consider an isolated system and calculate the residual energy that remains in the system at the end of a switching event.This is the energy that will need to be eventually dissipated.For quantum switching, there is a remarkable exponential reduction 45 in the residual energy with switching time: E ∝ e −aT s .For a rigid molecule, this has been extended to calculate dynamically the flow of energy between the molecule and the environment.For rapid switching, the power dissipation shows the characteristic quantum exponential decrease as the switching time is increased.For intermediate switching speeds there is an inverted region in which power dissipation actually increases with slower switching speed because of excitation from the thermal environment.At large switching times the power dissipation follows the classical 1/T s dependence. 46ere we extend the model to include the effect of atomic reorganization.Section II describes the mathematical model for the system and interaction with the environment.Section III shows the impact of reorganization energy on polarization and switching dynamics.We solve the switching dynamically using the Lindblad equation and observe hysteresis in cell polarization when the bias sweep direction changes.Thus we demonstrate that a single mixed-valence molecule can function as a memory element.In Section IV we calculate the power dissipated by switching a bit.A higher reorganization energy causes an increase in dissipated power.

II. TWO-STATE MODEL
A two-state quantum-dot system is illustrated schematically in Figure 3.It is the simplest QCA element consisting of two dots and a mobile charge.Charge transport between dots is through quantum-mechanical tunneling.We use |L and |R as the basis states and the binary digital data in the two-state QCA is encoded in the polarization P. When the charge is fully localized on the left dot, the system is in state |L and P = 1 representing a binary 0. When the charge is fully localized on the right dot, the system is in state |R and P = −1 representing a binary 1. FIG. 4. Schematic representation of charge transfer during switching.At t = 0 system is in ground state and charge is localized on left dot L. Bias ∆(t) is applied to the left dot and is linearly increased from ∆ initial to ∆ f inal over the switching time T s .Midway through the switching event at crossover t = T s /2, the energies of the two dots are matched.At the end of the switching event t = T s , charge is localized on the right dot and the electron transfer is complete.
Charge is transferred by changing the relative bias between the two dots.A simple switching operation is shown in Fig- ure 4. At t = 0, the charge is localized on the left dot.A time varying bias ∆(t) is applied to the left dot and the right dot is held at 0. ∆(t) is linearly increased from ∆ initial to ∆ f inal over switching time T s .Midway through the switching event at crossover t = T s /2, the energies of the two dots are matched.At the end of the switching event t = T s if the bias ∆ is large enough, the charge is almost entirely localized on the right dot and the electron transfer is complete.Polarization changes continuously from a polarization P = 1 at t = 0 to a polarization of P = −1 at t = T s in this example.We interpret polarization near 1 as a binary 0 and polarization near -1 as a binary 1.The density operator ρ can be written in terms of Pauli operators as where σi = Tr ρ σi .
The density operator is Hermitian with unit trace and the vector σx , σy , σz will lie on or within a unit sphere (commonly called the Bloch sphere): The Hamiltonian operator for the simplified electronic twostate subsystem is written as The characteristic time scale for dynamics is The natural scales for the problem are to measure energy in units of γ and time in units of T γ .
The polarization P is the expectation value of σz When the electronic charge occupies a dot, the atoms around the dot shift their positions in response.This lowers the total energy for the electron being on that dot.The energy associated with this relaxation of nuclear positions is the reorganization energy of Marcus theory λ . 47,48The relaxation is modeled as the lowering of the energy of the dot when occupied by the electron (as in Figure 2).
The Hamiltonian for linear coupling 49 between the electron and the motion of the surrounding ligands can be written in terms of λ : The anti-symmetric oscillating motion of the left and right dot ligands we treat as a classical harmonic driven by the imbalance between the electron occupancy between the left and right dots.The Hamiltonian operator representing this energy stored in the displaced ligand positions (classically 1 2 kx 2 ) is written as Accounting in this way for energy as it moves back and forth between the electronic and vibrational degrees of freedom guarantees conservation of total energy.The Hamiltonian is the sum of the Hamiltonians of the electronic subsystem ĤE , the electron-ligand interaction ĤEL , and the ligand subsystem ĤL The Hamiltonian depends on the electron state through the expectation value σz .This non-linearity is fundamentally due to the fact that the effect of the ligand distortion on the electron system has been treated semiclassically.
When a system is in equilibrium with an environment at temperature T the steady state density operator ρss can be written as The steady state Hamiltonian Ĥss = Ĥ( ρss ) depends on the steady state density operator ρss through Eq. (10).For given values of γ, ∆ (which may vary in time), λ , and T , Eq. ( 11) must be solved self-consistently.We will see that in some circumstances multiple solutions are possible for the same parameters.
For an isolated system, the time evolution of the density operator ρ is unitary and can be solved using quantum Liouville equation: An open system interacts with the environment and the time evolution of the density operator is non-unitary.For Markovian systems, under reasonable assumptions, the density operator can be solved using a Lindblad equation: [50][51][52] where { Â, B} is the anti-commutator of operators Â and B. The first term in Eq. ( 13) governs the unitary evolution of the density operator and is identical to Eq. ( 12).The Lindblad dissipator D models the system-environment interaction including dissipation and decoherence.The Lk are the Lindblad operators.We model the coupling of the system to the thermal environment for our two-state system as: Here Ĥ(t) |u k (t) = E k (t) |u k (t) , k B is the Boltzmann constant, and T is the temperature of the environment (bath).T d is a relaxation time or dissipation time and is a measure of the coupling strength between the system and the environment.T d is small for a system coupled strongly to the thermal bath.For an isolated system with no interaction with the environment T d = ∞ and both L1 and L2 are zero.The Lindblad operator L1 in Eq. ( 14) models the decay of the system from the instantaneous excited state to the instantaneous ground state by dissipating energy to the environment.The Lindblad operator L2 in Eq. ( 15) models the thermal excitation of the system from the ground state to the excited state by absorbing thermal energy from the environment.The operators are constructed such that when the potential driver stops at t = T s , the density operator will relax with characteristic time T d to the thermal equilibrium density operator in Eq. (11).

III. POLARIZATION OF AN OPEN TWO-STATE SYSTEM DURING SWITCHING A. Hysteresis in Polarization during Dynamic Switching
We calculate the polarization during a switching operation of an open two-state system with a Hamiltonian defined in Eq. (10).We solve the Lindblad equation in Eq. ( 13) for the density operator ρ(t) as a function of time and the polarization P is calculated directly from the density operator using Eq. ( 6).The system is in contact with the environment at temperature T and the coupling between the system and environment is characterized by the dissipation time T d .Initially at t = 0, the bias is ∆ min = −25γ and the system is in equilibrium with the environment.The initial dot occupancy probabilities are determined by the steady state density density operator from Eq. ( 11).The bias is linearly increased from ∆ min = −25γ at t = 0 to ∆ max = +25γ over switching time t = T s which we define as a positive sweep.After the positive sweep, the bias is held at ∆ max for a time sufficiently long such that system can come to equilibrium with the environment.Then the bias is linearly decreased from ∆ max back down to ∆ min over the same switching time T s which we define as a negative sweep.We explore the effect of changing the reorganization energy λ , switching time T s , dissipation time T d and environmental temperature T .Figure 5 shows the effect of the reorganization energy λ for quasi-adiabatic switching (T s /T γ ≫ 1).When there is no reorganization energy (λ = 0), the polarization of the positive sweep and the negative sweep overlap.At zero bias (∆=0), the polarization is zero and therefore no binary information is stored in the two-state system.
When the reorganization energy is introduced, the system exhibits hysteresis.The polarization is different depending on the sweep direction.The hysteresis widens with increasing reorganization energy.At t = T s /2 when the bias is zero (∆ = 0), the polarization P ≈ 1 or P ≈ −1 depending on the sweep direction.This means the binary information ("0" bit or "1" bit) is stored even when the bias is removed.The polarization is closer to 1 or -1 when the reorganization energy increases in-dicating bit storage is easier at higher reorganization energies.We will discuss this more in the next subsection.The slope of the polarization becomes more abrupt due to the effect of the reorganization energy.The onsite energies of the left dot E L and right dot E R are defined as In Figure 6, we plot the difference of the onsite energies of the left and right dots (E L − E R ) for the quasi-adiabatic switching calculations shown in Figure 5.When λ = 0, E L − E R simply matches the applied bias.The difference linearly increases from ∆ min to ∆ max during the positive sweep and linearly decreases from ∆ max to ∆ min during the negative sweep.Now consider the system with a nonzero reorganization energy.Initially when charge is localized on the left dot (P = 1) the minimum value of E L −E R is ∆ min −λ .During the positive sweep as the bias increases, E L −E R increases at the same rate.When E L comes closer to E R , the charge transfer begins.Increasing charge occupation causes relaxation of the nuclei in the right dot and lowers the energy of right dot.The lowering of energy accelerates the charge transfer which further lowers the energy of the right dot.This phenomenon causes switching to be less adiabatic and more abrupt.This is evident from the steeper slope of E L − E R in a system with reorganization energy compared to a system with λ = 0.The polarization curves in Figure 5 become abrupt with increasing reorganization energy.At the end of the positive sweep, the charge is localized on the right dot (P = −1) and E L − E R is ∆ max + λ .The same behavior is observed during the negative sweep but in the opposite direction.
Figure 7 shows the impact of switching time T s .When the switching is slow (T s /T γ = 1000) and therefore quasiadiabatic switching, the system follows the instantaneous ground state in accordance with adiabatic theorem.The charge switching is complete and the polarization changes from -1 to 1 and vice versa.When the switching is rapid (T s /T γ is small), the charge transfer is not complete at end of the sweep (t = T s ) and we observe charge quantum oscillations between the two dots post-crossover.The system will return to ground state a sufficiently long time after the switching event.
Figure 8 shows the impact of coupling strength between the system and the environment.When the system is strongly coupled to the environment (T d is small), any excess energy is immediately dissipated and the switching is more abrupt.For systems that are weakly coupled to the environment (T d is large), the system takes longer to reach the ground state and the complete switching of the charge takes longer.
Figure 9 shows the impact of environmental temperature.At higher environmental temperature, the hysteresis decreases and the switching becomes less abrupt.The system takes longer to reach the ground state due to thermal excitation from the environment.
From Figures 5 and 9, we observe that the width of the hysteresis is set by the reorganization energy λ and environmental temperature T .Changing either the switching time T s or the dissipation time T d (in Figures 7 and 8), does not change the hysteresis but only changes the time taken to reach steady state.

B. Memory Operations on a Two-state system with Reorganization energy
The hysteresis observed in Figure 5 shows that a system with nonzero reorganization energy can store a bit ("1" or "0") when the bias is zero during switching.We now perform a memory operation where we write and hold a "1" bit and then write and hold a "0" bit.The waveform of the applied bias is shown in Figure 10a and the polarization response is shown in Figure 10b.We split this operation into four regions (I to IV) for the ease of explanation.Consider the polarization response when λ = 5γ.FIG. 10.Writing and holding 1 and 0 bits for different reorganization energy λ .In a system with sufficient reorganization energy, the polarization does not change significantly when the bias is removed.The system continues to store the bit (1 or 0) that is written.When λ = 0, the polarization becomes zero when the bias is removed and the bit is not stored.Here λ /γ = [0, 5, 10], T s = 1000T γ , T d = 10T γ and k B T = 1γ.
In the region labeled I, the "1" bit is written.Initially at t = 0, the applied bias ∆ = 0 and the system is in equilibrium.The bias linearly increases and charge starts transferring.The bias is increased to ∆ = +25γ by which point the charge is completely localized on the right dot with polarization P = −1.The system now represents a binary "1."The bias is then held constant for some time so the system can settle to a equilibrium.The bias is now linearly decreased from ∆ = +25γ to ∆ = 0.In the system with significant reorganization energy (λ = 5γ), the charge does not transfer back and the polarization continues to be approximately -1.The binary information written into the system is locked.The bit is stored even when the bias is removed because the nuclear relaxation has localized the charge on the right dot.
In region labeled II, the bias is removed for a time duration of 10T d .The polarization for λ = 5γ does not change significantly and the binary information is not lost.The polarization continues to be approximately -1.
In the region labeled III, we write and hold a "0" bit and therefore the applied bias is of opposite sign of the waveform in region I.The bias is decreased from ∆ = 0 to ∆ = −25γ which causes the charge to localize on the left dot (P=1) and the system now represents a binary "0".When the bias is increased from −25γ to back to 0, the polarization does not change significantly for λ = 5γ.The polarization continues to be approximately equal to 1.
In the region labeled IV, the bias is held at 0 for a time duration of 10T d similar to region II.When λ = 5γ, the polarization continues to be approximately equal to 1 and the 0 bit is stored.
Figure 10 also shows the impact of reorganization energy.When λ = 0, there is no nuclear relaxation and the polarization is zero (bit not stored) in the absence of any applied bias.It is easier to hold binary information when the system has higher reorganization energy.

C. Multiple Steady State Configurations
In III A, we calculated the polarization of a switching system dynamically using the Lindblad equation.We observed different values of polarization for the same value of applied bias.This means the system can be in more than one stable states.We now calculate the steady state solutions of the system using Eq. ( 11) to determine under what circumstances the system will support a single or multiple equilibrium solutions.
We calculate the steady state density operator ρss at a given bias by self-consistently solving Eq. ( 11) using the constraint in Eq. ( 3).Without the loss of generality we can set σy = 0. To solve Eq. ( 11) iteratively, we start with an initial guess for the density operator and evaluate the RHS of Eq. ( 11) to obtain a value for ρss .This value is then substituted back into (11) and another iteration is obtained.These iterations continue until we obtain a converged solution for the density operator that satisfies Eq. (11).We then calculate the polarization from the steady state density operator using Eq. ( 6).
We calculate the steady state density operators for different values of bias ∆.Figures 11a-c shows the results for bias ∆ = 10γ, 2γ, 0 respectively.Here we take reorganization energy λ = 5γ and environmental temperature k B T = 0.25γ.The initial guess for the density operator is randomly chosen to lie on or within the unit circle defined by σx and σz .We then solve Eq. ( 11) iteratively to obtain a converged steady state solution for density operator and σx and σz .We perform this calculation for 100 random initial guesses for each bias condition by solving Eq. ( 11) for each of the intial guesses.In Figure 11a, b, c  Figure 11a shows the results when the applied bias is large (∆ = 10γ).After solving Eq. ( 11), all the initial random states yield a single solution for ρss with a polarization σz = −0.99 and the charge is localized on the right dot.
Figure 11b shows the results when the bias is lowered to ∆ = 2γ.All the initial random states converge to two possible solutions with polarization σz = −0.96and 0.73.The two polarization solutions are of opposite signs but not symmetric across zero.At a low applied bias, the system is bistable.For these steady state solutions, the majority of the charge can be localized on either dot.
Figure 11c shows the results when ∆ = 0.There are three possible steady states solutions (labeled A, B and C).Solutions A and B with σz = −0.92 and σz = +0.92 has the charge localized on the right or left dot respectively.Solution C with σz = 0 has the charge equally distributed between the dots.At a bias value close to 0, there can be three steady-state solutions.
Figure 11d shows the steady state polarization as a function of bias ∆.Multi-state equilibrium solutions are observed when the bias is comparable to λ + k B T .The polarization curve exhibits hysteresis.There are solutions near the origin (example solution C) that are not observed in the hysteresis behavior during dynamic switching calculations in Figure 5.To understand why, we calculate the expectation value of energy for solutions A, B and C.
In Figure 11e we calculate the expectation value of energy E for the three solutions A, B, C at ∆ = 0. E for the two symmetric solutions A and B are equal and E = −1.45γ.For solution C, σz = 0 and the expected value of energy E = −0.99γ.Solution C is not a stable state due to higher E compared to solutions A and B. Solution C will not be visited during dynamic switching.The system will select either solution A or B depending on the sweep direction.

IV. ENERGY DISSIPATION DUE TO SWITCHING A. Isolated System
We can first calculate the energy dissipated by considering an isolated system, uncoupled to the environment, and calculate the excess energy E excess left in the system at the end of the switching event (t = T s ).We know that for an open sys-tem this is the energy that will be eventually dissipated.What this approach ignores is dissipation, or thermal excitation, that occurs to the environment during the switching process.For all calculations reported in this section, ∆ is linearly increased −25γ to +25γ over switching time T s .This is sufficient to strongly localize the charge both before and after switching.
Consider the switching event described in Figure 4. We define the excess energy of the system E excess as the difference between the expectation value of the energy E and the ground state energy E 1 when the potential driver has stopped changing at t = T s .
The dimensionless adiabaticity parameter β captures the effect of parameters that impact the adiabaticity of the switching.
For an isolated system with no reorganization energy λ = 0, the excess energy E excess decreases exponentially with β There is no fundamental lower limit to E excess as it can be decreased to any arbitrary low value by increasing the switching adiabaticity β .The adiabaticity β can be made larger by switching more slowly (larger T s ), or using a smaller potential driver (smaller ∆), or by increasing γ.
The excess energy E excess left in the system at the end of the switching event is calculated using Eq. ( 18) for different reorganization energies λ and plotted as a function of switching time T s in Figure 12.In the absence of reorganization energy λ = 0, E excess decreases exponentially with increasing T s (and therefore β ) as shown in Eq. (20).This result has been discussed previously. 45The excess energy increases with the introduction of a nonzero reorganization energy λ .The reorganization energy causes the switching to be less adiabatic (see explanation for Figure 6 in the previous section).For the same switching speed, a system with higher reorganization energy has higher excess energy.

B. System in contact with a thermal environment
We now calculate E diss , the total energy dissipated for a system in contact with a thermal environment due to a switching event.We first calculate the energy dissipated E switch during switching from t = 0 to t = T s .We do this by calculating the power flow and then integrating it over t = 0 to t = T s .
The total power flow P total into the system is the derivative of expected energy of the system: FIG.12. Excess energy E excess at the end of a switching event for an isolated system as a function of the switching time T s for different values of the reorganization energy λ .The inset is the same data with the horizontal axis plotted on a linear scale.E excess decreases exponentially with T s when λ = 0. Increasing λ causes switching to be somewhat less adiabatic and increases E excess .
Using the Lindblad equation in Eq. ( 13), we can write The first term vanishes because of the cyclic property of the trace and Eq. ( 24) reduces to By expanding the second term in the above equation using Eq. ( 9) we obtain (26) By substituting Eq. ( 7) in the third term in Eq. ( 26) we get So Similarly using Eq. ( 8) we obtain The third term P 3 and fourth term P 4 in Eq. ( 26) cancel each other and therefore the total power flow into the system reduces to We identify the second term P 2 in Eq. ( 30) as the power flowing into the system due to work done by the time-varying control electrodes during switching.
We identify the first term P 1 in Eq. ( 30) as the instantaneous power flow into the system from the environment.The power flowing from the system into the environment is obtained by simply changing the sign.
The total energy dissipated during the switching event E switch is the integral of instantaneous power dissipated P switch during entire switching event from t = 0 to T s .
We now consider the residual excess energy E excess left in the system when the potential driver stops at t = T s .This residual excess energy will be dissipated after some time t > T s .We define E excess as the difference between the non-equilibrium expectation value of the energy and the steady-state expectation value at the end of the switching event. where E ss (T s ) = Tr ρss (T s ) Ĥss (T s ) The steady state values of ρss and Ĥss are obtained by solving Eq. ( 11) self consistently.The total energy dissipated E diss is the sum of the energy dissipated during switching E switch and excess energy E excess left at the end of switching event which is eventually dissipated to the environment.E diss = E switch + E excess (37)   Substituting Eqs.(33), (34), (36) and (36) in the above equation, we obtain For an open system we solve for the density operator ρ using the Lindblad equation in Eq. (13).We then calculate the total dissipated energy E diss using Eq. ( 38).E diss is calculated as a function of switching time T s for different reorganization energies λ .We first focus on the λ = 0 case shown in Figure 13, which is our previously reported result. 46The E diss vs. T s behavior can be broken into three regimes.
First, the exponentially adiabatic regime at rapid switching (low T s ) which is labeled I.When the switching is rapid, there is significant excess energy E excess left in the system at the end of the switching event and the amount of energy dissipated during switching E switch is negligible.The total dissipated energy E diss ≈ E excess and decreases exponentially with switching time as shown in Eq. (20).
Second, the inverted regime at intermediate switching speeds labeled II.In this regime, we observe the total dissipated E diss increases as the switching speed decreases.It is surprising that slower switching produces higher dissipation.At these intermediate switching speeds, E excess is negligible due to its exponential dependence on T s .The total dissipated energy is approximately the energy dissipated during switching E diss ≈ E switch .When the switching time is much faster than the dissipation time constant T d , E diss is small because the switching is complete before any thermal excitation from the environment can occur.As the switching speed slows down, thermal fluctuations from the environment cause excitations above the ground state.In this regime, the system is moving slowly enough for the excitations to occur, but not slowly enough for them to immediately de-excite and permit the system to closely track the equilibrium occupations.The result is that slowing the switching speed counter-intuitively causes E diss to increase.We have shown that the inverted region does not occur when environmental temperature is zero.The peak of the inverted region occurs at a switching time T s ≈ 2T d .This inverted dependence has been confirmed using a semiclassical treatment. 46hird, the classical regime which is labeled III, at slow switching speeds (large T s ).Here the total dissipated energy E diss decreases linearly with T s matching a classical system.At slow switching, the system is close to equilibrium with the environment and any thermal excitation or de-excitation occurs rapidly.But the system is slightly excited and this energy decreases linearly with increasing switching time.The dotted line in Figure 13 showing the 1/T s dependence has been added as a visual guide.
The impact of reorganization energy is shown in Figure 14.When the reorganization energy is small (λ /γ 1), we observe exponential dependence at rapid switching (small T s ).But the rate of E diss decrease is slower.At rapid switching the dissipated energy is approximately equal to the excess energy left at the end of switching event E diss ≈ E excess .This excess energy increases with increasing reorganization energy as shown in Figure 12.As switching slows down, E diss matches system with no reorganization energy (λ = 0) and we observe the inverted region and the classical dependence.
At intermediate reorganization energies (1 λ /γ 3), the dissipated energy increases in the exponentially adiabatic region.At slower switching E diss meets the inverted dependence at λ = 0.The inverted region gets less pronounced with increasing reorganization energy.
At large reorganization energies (λ /γ 3), the dissipated energy is very high at rapid switching.At slower switching, the inverted region completely vanishes and the dissipated energy follows the classical 1/T s dependence.
Figures 15 and 16 show the impact of environmental coupling when reorganization energy λ = 1γ and λ = 6γ respectively.T d /T γ is varied by four orders of magnitude from 10 −1 to 10 3 at environmental temperature k B T = 3γ.The dotted line shows the 1/T s dependence of a classical system.When the system is very strongly coupled to the environment (T d is small) any excess energy generated during switching is immediately dissipated and the system acts classically.The system does not exhibit the characteristic quantum adiabatic exponential decrease even at fast switching which is evident in the T d /T γ = 0.1 case which shows E diss linearly decreasing with T s .For small reorganization energy λ = 1γ as the coupling becomes weaker, E diss follows the exponential decrease  13.For small values of λ , E diss is higher in the exponential adiabatic regime.E diss is similar to E diss for λ = 0 in the inverted and classical regimes.For larger values of λ , E diss is increased in the exponential adiabatic region compared with the λ = 0 case and the inverted region vanishes.At higher T s , E diss shows a classical dependence on T s .Here the bias ∆ is varied from −25γ to +25γ over T s , T d = 100T γ and kT = 3γ.The dotted line is the 1/T s dependence of a classical system. of the isolated system at lower T s , but eventually, as T s increases, the interaction moves dissipation into the classical viscous regime.The inverted region between the two regimes is readily apparent.For larger reorganization energy λ = 6γ, E diss decreases linearly with T s and increases as environmental coupling weakens.
Figure 17 shows the results for different environmental temperatures with T d = 100T γ .We observe that when switching is rapid, E diss matches the excess the energy of an isolated system.When the temperature is low there is no significant thermal excitation from the environment and E diss decreases monotonically with T s .As the temperature increases E diss starts exhibiting the inverted regime at larger switching times.

V. DISCUSSION
Molecular relaxation and subsequent lowering of dot energy upon charge occupation has been reported in candidate QCA molecules.We have here examined the impact of this molecular reorganization energy on the localization and energy dissipation of a two-state QCA cell.
The Hamiltonian of a QCA cell with reorganization energy, Eq. ( 10), is non-linear due to the effect of ligand distortion.The steady state Hamiltonian depends on the steady state density operator and for small applied biases can have multiple solutions for the same parameters, as seen in Figure 11.We solved for the switching dynamics using the Lindblad equation (13) to model the environmental interaction at finite temperature.Molecular QCA cells with non-zero reorganization energy display hysteresis in the polarization response, as shown in Figure 9.
A stored memory bit must necessarily be in a long-lived metastable state that depends on the system's past.Because molecules cannot be addressed individually, the usual QCA approach is to hold binary information not in a single cell, but rather in moving bit packets comprising several cells that are shepherded smoothly around a circuit by the clocking fields. 22,53The device feature size can be much smaller than the electrodes that provide the clocking field which move information in computational waves around the circuit architecture. 24,54The stability of the bit is then due to the large kinetic barrier to flipping all the cells in a bit packet, a stability enhanced by the quantum decoherence caused by the environment. 55he reorganization energy λ stabilizes the localized charge when the bias is removed, allowing a single molecular twostate QCA cell to act as a memory cell for storing a single bit.The molecular reorganization energy provides another source of bit stability, namely the kinetic barrier associated with rearranging the atoms within the molecule in addition to moving the electron from site to site.The reorganization energy can be larger or smaller than k B T , depending on the details of the molecule and its connection to the surrounding environment.
The enhanced bit stability due to the intrinsic molecular reorganization energy comes at the expense of somewhat higher power dissipation for a switching operation.As we have seen in Figure 6, the movement of the dot ligands in response to the electron transfer accelerates the electron transfer process, in a way reminiscent of feedback.The result is lower adiabaticity and more energy dissipation for the same switching speed.It is to be noted that this effect is more pronounced at slower switching speeds than at the higher speeds.If the switching speed is faster than the speed at which the ligands can respond, then the reorganization is irrelevant.For small reorganization energies, we observe three distinct regimes in dissipated energy dependence on switching time: an exponential adiabatic (e −aT s ) regime, an inverted regime in which dissipation increases with slower switching, and a classical 1/T s dependence for very slow switching.For large reorganization energy, the power dissipation is primarily classical.
From a practical molecular QCA design perspective, one has several key parameters to consider.The tunneling en-ergy γ and therefore T γ is determined by the choice of the linker between the charge centers.The coupling strength to the environment T d depends on the precise way the molecule is bonded to substrate-essentially the amount of heat-sinking.The switching signal E c applied through electrodes and switching time T s can be modified through the external driver circuitry.A mixed-valence molecule in a polar solvent typically has a large reorganization energy because many surrounding molecules must move to respond to an electron transfer event.For a surface-bound molecule not in solution, this energy can be much less.We have explored a range of λ in our calculations.
The results show a trade-off between enhanced QCA cell bi-stability and energy dissipation.One might expect that strong coupling to the environment, i.e., strong heat sinking (low T d ), would be optimal, but our results show otherwise.Systems with modest coupling to the substrate can take advantage of exponential adiabaticity at fast switching speeds (small T s ) and therefore optimize net power dissipation.
DATA AVAILABILITY

FIG. 1 .
FIG. 1.Chemical structure of the mixed-valence diferrocenylacetylene (DFA) molecule.The iron centers in the two ferrocenyl groups act as dots for charge localization, forming a molecular double-dot.An extra electron present on one or the other dot encodes binary information.The energy associated with the relaxation of each ferrocenyl group in response to the presence or absence of the electron is the reorganization energy.

FIG. 5 .
FIG. 5. Hysteresis in polarization and the impact of reorganization energy on an open system.Increasing reorganization energy λ increases the hysteresis in polarization and causes switching to be less gradual.With sufficient λ , polarization is approximately equal to +1 or -1 at zero bias (t = T s /2) depending on the bias sweep direction.In this calculation, the bias ∆ linearly increases from −25γ to +25γ over switching time T s during the positive sweep and decreases from +25γ to −25γ over the same switching time T s during the negative sweep.Here λ /γ = [0, 5, 10], T s = 1000T γ , T d = 10T γ and k B T = 1γ.

FIG. 6 .
FIG. 6. Difference in onsite energy of dots E L − E R for the switching operation shown in Figure 5.With no reorganization energy λ = 0, E L − E R matches the applied bias ∆.With a non-zero reorganization energy, the change in E L − E R is faster (steeper slope) compared to the change in the applied bias.This causes the switching to be more abrupt and therefore less adiabatic.

FIG. 7 .
FIG. 7. Hysteresis in polarization and the impact of switching time for an open system.For quasi-adiabatic switching (large T s /T γ ), charge transfer is smooth and complete.Switching too rapidly causes quantum oscillations and incomplete charge transfer.The system may be left in excited state at the end of the switching event.In this calculation, the bias ∆ linearly increases from −25γ to +25γ over switching time T s during the positive sweep and decreases from +25γ to −25γ over the same switching time T s during the negative sweep.Here T s /T γ = [10, 100, 1000], λ = 10γ, T d = 10T γ and k B T = 1γ.

FIG. 8 .
FIG. 8. Hysteresis in the polarization and the impact of coupling strength to the environment for an open system.Systems coupled strongly to environment (small T d ) reach steady state quickly.Weakly coupled systems (large T d ) take longer to reach steady state.In this calculation, the bias ∆ linearly increases from −25γ to +25γ over switching time T s during the positive sweep and decreases from +25γ to −25γ over the same switching time T s during the negative sweep.Here T d /T γ = [0.1, 10, 100], T s = 1000T γ , λ = 10γ and k B T = 1γ.
the random initial guesses are shown using open circles and the converged solutions for each of these initial guesses are shown using filled circles.

FIG. 11 .
FIG. 11.Multiple equilibrium solutions of the density matrix for different applied biases.Solutions for the steady-state density operator ρss are obtained by randomly choosing an initial state and solving Eq. (11) self-consistently for λ = 5γ and k B T = 0.25γ.In (a) to (c), the initial guesses are shown using open circle markers and steady-state solutions are show using filled circles.Steady-state solutions are calculated for three different bias values: (a) ∆ = 10γ.All initial guesses converge to a single steady-state solution.(b) ∆ = 2γ.There are two steady-state solutions and the system is bistable.(c) ∆ = 0.There are three steady-state solutions which are labeled A, B and C. (d) Steady-state polarization as a function of the applied bias.For small bias values there are multiple steady-state solutions.(e) Expectation value of energy E for the three solutions A, B, and C at ∆ = 0 shown in (c) and (d).E for solutions A and B is lower than E for solution C. The result is that solution C is not visited during dynamic switching.

FIG. 13 .
FIG. 13.Total dissipated energy E diss of an open system as a function of switching time T s with no reorganization energy (λ =0).For rapid switching (small T s ) in the regime labeled I, E diss exhibits exponential adiabaticity.For intermediate values of T s in the regime labeled II, we observe the inverted region where E diss actually increases as switching slows down.For large T s in regime labeled III, E diss follows the classical 1/T s dependence.Here the bias ∆ is varied from −25γ to +25γ over T s , T d = 100T γ and kT = 3γ.The dotted line simply indicates the slope of a 1/T s dependence on switching time that is characteristic of a classical system.

2 FIG. 14 .
FIG.14.Total dissipated energy E diss of an open system as a function of switching time T s for different reorganization energies λ .The λ = 0 calculation is the same as shown in Figure13.For small values of λ , E diss is higher in the exponential adiabatic regime.E diss is similar to E diss for λ = 0 in the inverted and classical regimes.For larger values of λ , E diss is increased in the exponential adiabatic region compared with the λ = 0 case and the inverted region vanishes.At higher T s , E diss shows a classical dependence on T s .Here the bias ∆ is varied from −25γ to +25γ over T s , T d = 100T γ and kT = 3γ.The dotted line is the 1/T s dependence of a classical system.

2 FIG. 15 .FIG. 16 .
FIG.15.Total dissipated energy E diss of an open system as a function of switching time T s for different values of the dissipation time constant T d when the reorganization energy λ is small.Here the bias ∆ is varied from −25γ to +25γ over T s , λ = 1γ.and kT = 3γ.The dotted line is the 1/T s dependence of a classical system.

2 FIG. 17 .
FIG. 17.Total dissipated energy E diss of an open system as a function of switching time T s for different values of the environmental temperature T .E diss increases as T is increased.The inverted region becomes more pronounced with increasing T.Here the bias ∆ is varied from −25γ to +25γ over T s , λ = 2γ and T d = 100T γ .