Exploring the write-in process in molecular quantum cellular automata: a combined modeling and first-principle approach

The molecular quantum cellular automata paradigm (m-QCA) offers a promising alternative framework to current CMOS implementations. A crucial aspect for implementing this technology concerns the construction of a device which effectively controls intramolecular charge-transfer processes. Tentative experimental implementations have been developed in which a voltage drop is created generating the forces that drive a molecule into a logic state. However, important factors such as the electric field profile, its possible time-dependency and the influence of temperature in the overall success of charge-transfer are relevant issues to be considered in the design of a reliable device. In this work, we theoretically study the role played by these processes in the overall intramolecular charge-transfer process. We have used a Landau–Zener (LZ) model, where different time-dependent electric field profiles have been simulated. The results have been further corroborated employing density functional tight-binding method. The role played by the nuclear motions in the electron-transfer process has been investigated beyond the Born-Oppenheimer approximation by computing the effect of the external electric field in the behavior of the potential energy surface. Hence, we demonstrate that the intramolecular charge-transfer process is a direct consequence of the coherent LZ nonadiabatic tunneling and the hybridization of the diabatic vibronic states which effectively reduces the trapping of the itinerant electron at the donor group.

S Supplementary material for this article is available online (Some figures may appear in colour only in the online journal) Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
Quantum cellular automata (QCA) have been proposed as an alternative and revolutionary paradigm to conventional semiconductor-based information processing architectures, not only for representing classical binary information using quantum mechanical systems but also to perform classical computation [1]. In the original paradigm, the basic computational cell consists of interconnected quantum dots so as to allow charge localization and charge exchange between pairs of dots [2]. Thus, the charges can move between strongly localized electronic states while controlled by an external electric field [3]. However, in this implementation this behavior is only preserved at low temperatures making it therefore impractical for actual computation. Attempts have been made to find alternative systems in which the paradigm can be successfully implemented [4]. In this sense, the molecular m-QCA concept has attracted a great deal of the attention offering the possibility to achieve room temperature functionality by developing tailored molecular complexes [5][6][7][8]. In the m-QCA implementation one of the crucial factors is the capability to effectively control the initial polarization state of the first molecular complex via an external electric field. Initially, an STM device has been employed simultaneously as a source of the external electric field and as a monitoring system to investigate the charge-transfer process [5]. The performed measurements demonstrated that the charge-transfer process causes an enhancement in the measured capacitance, which can be uniquely ascribed to the creation of an instantaneous molecular dipole moment that can be differentiated from the initial situation. Recently, a second experimental setup has been proposed in which the first molecular complex is placed between a nanogap of gold electrodes, where a voltage-drop can be generated and used as a source for the external electric field [9]. In both implementations (schematically shown in figure 1 for the purpose of illustration) the creation of the external electric field can be regarded as a dynamical process in which the driving force is actively tuned creating a timedependent profile rather than a homogeneous constant field.
In this work, we develop a theoretical model in which we study the dynamical nature of the charge-transfer process as a function of a time-dependent external field. Firstly, a Landau-Zener (LZ) model is constructed to investigate the dynamical change of the electronic occupation number and the influence of the nonadiabatic corrections in the overall efficiency of this process [10,11]. Subsequently, we have extended the model to describe nonlinear field profiles to address the influence of the time-dependent profile on the charge-transfer process. Our results have been further validated for all considered time-dependent electric field profiles by employing an in-house version of the DFTB + code [12,13], in which the different external field profiles have been coupled to the time-step of a first-principle based molecular dynamics simulation. The developed methodology has been applied to an alkyl-derivative C 37 H 50 N 4 O 4 (hereafter referred to as PNX) in which two localized electronic states are found in the highestoccupied-molecular-orbital (HOMO) and lowest-unoccupiedmolecular-orbital (LUMO) molecular states [14,15]. Finally, we address the influence of the time-dependent external electric field (within the linear sweep rate approximation) in the nuclear motion of the molecule by solving the complete Schrödinger equation where nuclear quantum effects are included. We show that the response of the electronic system to the time-dependent field creates a change in the adiabatic potential energy surface (PES) generating a double-well bistable potential energy profile, a result which highlights the important role played by this effect as a complementary mechanism that enhances the efficiency of the charge-transfer process.

Electronic structure calculations
We consider the simplest basic m-QCA (half) cell represented by the PNX complex, see figure 2. This molecule is composed of a dinitrobenzene moiety, which is an electron acceptor group, and the dihydrophenazine sub-unit which is a donor electron complex. The two moieties are bridged by a methylene group allowing us to control the spatial separation between donor and acceptor functional groups [15]. The geometry optimization and the corresponding electronic structure calculations have been carried out using the DFTB + code [16][17][18][19] and contrasted with the results obtained by employing all-electron density functional theory (DFT) at the B3LYP/6-311 + +G(d,p) theory level as implemented in the Gaussian 09 code (the results are presented in the SI (stacks.iop.org/ JPhysCM/31/405502/mmedia)) [20,21]. Although the bandgap energy is underestimated the overall trends are preserved leading to a valid qualitatively correct results.
The HOMO (figures 2(a) and (d)) and LUMO (figures 2(b) and (e)) frontier molecular orbitals for the n = 3 and n = 5 bridge lengths in the PNX molecular complex clearly displayed the previously mentioned electronic localization features, in agreement with previous works [14] .The computed values of the dipole moment are 5.45 D in the case of n = 3 and 5.7 D in the case of n = 5 with the strongest component along the direction of the donor-acceptor axis. All quantities have been obtained using the DFTB methodology [16][17][18][19].
This can be easily realized in the charge distribution difference displayed in figures 2(c) and (f) where a depletion of electrons is observed in the Nitrogen atoms, while an increment of charge density is seen at the Oxygen atom belonging to the Nitrogen dioxide moieties. Vibrational analysis for both moieties has been carried out and no imaginary frequencies were obtained allowing us to conclude that the optimized geometries lie on local minima in the corresponding multi-dimensional PES. Finally, the suggested molecules can be regarded as an expected prototypical moiety in which the m-QCA can be implemented where two well-localized and spatially separated electronic states can be uniquely assigned to the '0' and '1' classical bits necessary to represent information.

A dynamical model for the switching behavior
To describe the dynamical switching behavior associated with the intramolecular charge-transfer process, the simplest approach relies on a two-state approximation where the timedependent Hamiltonian can be written as: Here, localized diabatic states are used to represent the initial and final stages of the charge-transfer process (the details of the approximation are explained elsewhere, see [22,23]). In our case, one of the charges is assumed to be transferred from the occupied HOMO to the unoccupied LUMO leading to the normalization condition: |c a | 2 + |c b | 2 = 1, where c a and c b represent the occupation number of the HOMO and LUMO molecular orbitals, respectively. The Hamiltonian matrix elements can be written as The corresponding eigenfunctions (in the diabatic representation) for H aa and H bb are denoted by |φ a and |φ b with associated eigenvalues E a and E b , respectively. The scheme of the two-state model is presented in the lower panel of figure 1. Following a similar line as in [23], we can recast the time-dependent diabatic states as: Here, L denotes the distance between driver and target and d is the intramolecular distance between the donor and acceptor moieties. At this point, two important quantities are introduced that describe the simultaneous dynamical change occurring in the the molecular complex (target) and the source of the external electric field (driver). Firstly, the time-dependent change in the dipole moment of the molecular complex can be followed by defining the function P 1 (t) = 2|c a (t)| 2 − 1 that indicates not only the change in the population of the involved energy levels of the target, but also defines a metric to characterize the switching behavior needed to compute. Similarly, one can express the change of the external electric field (driver) in terms of the reorganization of charge within the non-invasive electrodes and define P 1 (t) = q 1 (t) − q 2 (t) assuming that the two partial charges q 1 2. An external driver such as an electric field (static or time-dependent) acts on a molecular system (target). Within a quantumclassical approach, the (classical) driver consists of two partial charges q 1 and q 2 satisfying q 1 + q 2 = 1 corresponding to one mobile electron. The (quantum) target is symbolically represented by two centers (triangles) a and b between which charge transfer can take place (donor and acceptor). The corresponding real occupation probabilities are c 2 a,b , so that a target polarization can be defined as P 2 = 2c 2 a − 1. The intramolecular charge transfer in the target can be mapped onto classical '1' and '0' bits of information. The distance between the electric field source and the molecular complex is denoted by L and the distance between donor and acceptor is denoted by d.
and q 2 created in the electrodes correspond to a single mobile electron. Noteworthy to mention is that one can also define diabatic states by employing the matrix representation of the electric dipole moment operator: [24,25] (refer to SI for a detailed discussion). The main advantage of this representation is the possibility to include geometrical factors of the alignment between the electrodes and the molecular complex. Subsequently, in order to obtain the ground state of the system one has to dynamically transform the diabatic basis set to the adiabatic basis set by employing a unitary transformation at every time-step of the propagation [26][27][28]. Once this is done one can write the time-dependent Schrödinger equation (TDSE) as: where E ± are the adiabatic energies and Ψ(t) represents the adiabatic wave functions, The unitary transformation U(t) defines a new basis set, the adiabatic states, such that (Ψ + , is introduced, related to a dynamical mixing of the diabatic electronic states (nonadiabatic corrections). At this point, one can use a similar methodology to the Jacobi eigenvalue algorithm [29] and obtain a relationship for the angle θ of the unitary transformation matrix as: Hence, the ratio at which the θ parameter changes over time is analytically obtained as: with dEz(t) dt as the rate at which the electric field changes over time [26]. In this case, the studied γ parameters belong to the PNX molecule with n = 3 and n = 5 units in the bridge and can be calculated by employing Koopman's theorem. The obtained values for γ are 0.45 eV and 0.27 eV for n = 3 and n = 5, respectively. The results of the simulations with the minimal model are displayed in figure 3, where the response function P 2 (t) and the θ parameter are first obtained as a function of a linearly varying time-dependent external electric field. Our result indicates, as expected, a better switching behavior in the n = 5 than in the n = 3 moiety under the same initial geometrical conditions [28]. Henceforth, we will carry out the rest of the studies employing the PNX n = 5 molecular complex, while the results for the n = 3 case are shown in the SI. In terms of the nonadiabatic corrections one can observe a non-linear increment (see insets in figure 3) of this parameter, especially strong in the time interval where the two-level system is passing through an avoided crossing. This result indicates the important role played by these corrections in the intramolecular charge-transfer process. Conceptually, the LZ model [10,11] offers a powerful theoretical interpretation in which one can understand the behavior of the system in terms of the internal time scales of the molecular complex (τ int = γ ) and the external time scale (τ ext = E+−E− ) supplied by the driver. Thus, since the external electric field is driving the system through an anti-crossing regime with a sweep velocity ( dE dt ) the adiabatic crossing will be obtained only when the condition τ external >> τ internal is fulfilled, leading to the conditions needed to invoke the quantum-adiabatic theorem [26]. A consequence of this interplay between internal and external time scales is that the enhancement of the nonadiabatic corrections is a fingerprint of the adiabatic charge-transfer process.

Nonlinear time-dependent electric field profiles
In order to study a plausible scheme for controllable LZ tunneling in a molecule, one can propose the use of different electric field profiles that will provide different sweep velocities and therefore increase or decrease the probability to obtain an adiabatic charge-transfer process [27,30]. Hence, in this section a set of three non-linear pulses have been selected with the following time-dependence [31]: where ω is the corresponding frequency of the pulses and τ is the characteristic time of the pulse in the Gaussian pulse. The different profiles are shown in figure S3 of the SI where the shape and duration of the pulse have been displayed. The obtained results indicate that the choice of a pulse strongly affects the non-linear response of the molecular complexes. This is clearly observed in the case of the exponential pulse,  in which the switching behavior is diminished as shown in figure 4). Likewise, the exponential and cosine pulses exhibit a leading and trailing edges that pass twice through the resonance region where the adiabatic approximation breaks down and the system undergoes transitions between the ground and excited states via the LZ mechanism.
Hence, this effect can be understood in terms of the behavior of the θ parameter that accounts for the nonadiabaticity of the system. In both situations, the electronic nonadiabatic correction displays two peaks, one upwards and one downwards that are directly related with the process of HOMO-LUMO hybridization and the double switching pattern observed in the response function, see figure 4. Finally, the tanh pulse presents better results in terms of charge-transfer stability and can be suggested as a controllable LZ electric field pulse since a single switching behavior is found after the system has been driven through resonance ( figure 4). Hence, one can state that the quest of looking for a suitable candidate to implement the m-QCA should not only be reduced to find optimal chargetransfer parameters (γ), but also to consider the shape and sweep rate of the external electric field as an integral part of the write-in mechanism.

Validating the minimal model: first-principle calculations
To validate the results of the model, first-principle calcul ations were carried out in the selected PNX molecular complex. The dynamical nature of the electric field has been simulated by coupling the time-dependence of the external electric field with the time-step of quantum molecular dynamics as implemented in the DFTB + code [16,17,19]. Moreover, an in-house extended version of the mentioned code has been developed where many different electric field profiles such as the ones used in the L-Z model can be utilized [12,13]. All simulations have been performed at T = 30 K and the molecule has been previously thermally equilibrated to this temper ature during 1 ns. Subsequently, the external electric field has been turned on with different profiles: linear, exponential, cosine and tanh.
The results of these calculations for PNX with n = 5 are displayed in figure 5. Figure S3 shows the n = 3 results. One can clearly see that the switching behavior is obtained where P 2 (t) −0.9 for all different time-dependent external electric field profiles. In the specific case of the linear profile, one finds the corresponding S-shape behavior observed in the L-Z model around the resonance field. Concerning the exponential and cosine profiles, the charge transfer is achieved within a few femtoseconds (exponential) or for longer period of times (cosine) depending on the duration of the pulse and its periodicity as predicted by our model. Finally, the tanh profile offers an interesting case in which the switching is fast achieved and is constantly maintained along the molecular dynamics trajectory. In any case, the discrepancies in the efficiency computed by the minimal model and the DFTB + calculations can be ascribed to the fact that geometrical relaxation and charge reorganization of the molecular complex are considered in the first-principle calculation while those effects are not taken into account in the minimal model. One can conclude that the profile of the external electric field plays an important role in the response behavior, offering an efficient way to control the intramolecular charge-transfer process.

Quasi-static potential energy surface
In the previous models only the electronic terms of the total Schrödinger equation (i.e including the nuclear kinetic energy term) have been considered to be affected by the external electric field. However, recent theoretical works have indicated the importance of the electron-vibration coupling by producing considerable self-trapping effect which significantly influences key characteristics such as the cell-cell response [32]. Thus, in order to understand the role played by this interaction, we have extended our model to explicitly include the electronvibration coupling by solving the complete Schrödinger equation to obtain the corresponding vibrational wave functions. The model suggested here describes a substantially different situation than in previously studied theoretical works, where the potential energy in which the nuclei vibrate is created by the driven electron transfer process and not by the small displacements along the normal mode coordinates [33,34]. Therefore, this model is intended to study the nuclei behavior when subject to this instantaneous electronic potential created by the interaction of the charge-carrier and the time-dependent external electric field. The Hamiltonian that describes this situation can be written as: where the time dependent term V e ef ( r, t) refers to the timedependent external electric field that exclusively affects the electrons [34]. In order to solve the TDSE one can define a time-dependent basis set that provides the instantaneous eigenstates of the instantaneous total Hamiltonian and that can be written as: in which we have used the spatially localized crude-adiabatic (CA) diabatic states ψ CA D/A,Q0 ( r), where the itinerant electron is localized either at the donor or acceptor sites of the molecule [35]. This new molecular wave function is based on the electronic diabatic wave functions which have the same character at all values of the coordinate Q. The symbol Q 0 is used to emphasize that we optimized the atomic coordinates at every external electric field value ensuring that we are working with the ground state configuration. Once the atomic optimization is performed, we compute the new set of vibrational modes which are subsequently used to calculate the strenght of the electron-vibration interaction providing the value of this parameter to be used in our model [36]. Thus, the corresponding eigenvalue equation can be recast as: where the eigenvalues ε QS i (Q, t) are parametrically dependent on the normal mode Q [33]. Within our two-state model we may write the Hamiltonian as: with the matrix elements given by the following expressions: Q represents the chosen (antisymmetric mode) dimensionless normal mode, Q m is a displacement that locates the two harmonic potentials at different nuclear geometries, ω is the vibration frequency of the harmonic diabatic oscillators and E 0 represents the energy change in the charge-transfer process and calculated as defined in [24]. Thus, the total wave function is approximated as a product of the relevant electronic state ψ CA i,Q0 and a nuclear wave function χ nt (Q) and can be recast as: whence the corresponding Schrödinger equation can be expressed as: with, (16) In this model, the parameter λ(E z ) = 2 ωQ 2 m (E z ) is the electric field dependent reorganization energy and γ is the coupling between the electronic diabatic states (CA) [33,37]. The W ii (Q,t) term represents the potential created by the associated Harmonic diabatic oscillator and the time-dependent external electric field. In order to calculate numerically the reorganization energy, one can resort to a linear approach where a Taylor expansion in the molecular electronic levels is carried out [38]. In this case, the expansion will be formulated in terms of the vibrational coordinates Q i , where i denotes the mode vector allowing to write the energy as: where the displacement coordinate x i along the mode vector i, the vibrational coordinates X i and the cartesian atomic coordinates R i have been introduced. Generally, as a first approximation one can assume that the first order correction is the most important contribution, allowing to write the dimensionless electron-vibration coupling as: The electron-vibration couplings for the selected molecule have been calculated by using the aforementioned approximation and the actual values obtained employing the DFTB + code [16,17,19]. The results of these calculations are presented in figure S4 in the SI. Subsequently, once the electron-vibration couplings have been obtained, the solution of the vibrational Schrödinger equation has been obtained by employing the Fourier grid Hamiltonian method [26] with 256 grid points along the mass-weighted normal coordinate Q and the results are displayed in figure 6. The choice of this specific numerical method is motivated by the discreteness of our PES which is evaluated at given grid points and yielding directly the bound states eigenvalues and eigenfunctions at every studied value of the intensity of the external electric field.
We focus on the case E(t) = −0.05, 0.0, 0.05, i.e. before, during and after the resonance field. The ground state vibronic eigenvalue is localized at one of the most protruding minima of the adiabatic PES when E(t) = −0.05 and 0.05, namely, after or before the charge-transfer process has occurred. This fact can be ascribed to the vibrational trapping effect for the associated electronic wave function state. Concerning the resonance field case, one can observe the appearance of a bistability behavior in the adiabatic PES and a delocalization effect of the total vibrational wave function. This indicates that once the electric field pulse has achieved the resonance value where the charge-transfer process takes place, the associated adiabatic PES displays a double-well structure that hybridizes the Donor and Acceptor vibrational wave functions allowing the electron to be transferred. This issue can be understood in terms of the delocalization of total vibrational wave function over the molecular complex leading to an effective reduction of the electron trapping effect and therefore opening an instantaneous channel in which the diabatic vibrational modes associated to the donor and acceptor sites are degenerated. Thus, one can state that despite the fact that the normal mode that facilitates the electron transfer within the molecular complex was not the strongest one coupled at zero electric field (see SI for the actual values), the continuous change of the time-dependent electric field modifies the geometry reference point, creating a favorable condition to enhance the electron-vibration coupling for this mode and, hence, easing the charge-transfer process.
Finally, our findings highlight an intricate interplay between the external electric field, the charge transfer process, and its influence in the shape of the adiabatic PES. In this sense, we would like to emphasize the role played by the electron-vibration coupling and the vibronic levels corresponding to the energy region of avoiding crossing (where the nonadiabatic transitions is taking place) as the parameters which give a major contribution in the change of the PES shape.

Conclusions
An effective write-in strategy for an ideal m-QCA setup has been theoretically studied by combining minimal-model and DFT-based approaches. In both cases we have taken into account time-dependent external electric fields with different profiles and the results of both approaches are in good agreement. The current investigation suggests that different electric field profiles are related to different responses of the first m-QCA cell and may have deeper implications in the subsequent propagation of information and, hence, potentially alter the network functionality. Further investigations have been carried out in order to understand the role played by the vibrational modes of the molecule in the chargetransfer process. Thus, we have demonstrated that there is a correlation between the driven charge-transfer process and the enhancement of a specific electron-vibration coupling. An important conclusion is that within a class of potential m-QCA candidates, not only is the electronic coupling a crucial parameter, but also the electron-vibration coupling to specific modes easing the ET process. In any real-world implementation of this paradigm, the write-in system will have to account for these effects and therefore clarifying the origin of these issues is an important task to optimize the final device. There are obviously many open questions left, especially the influence of thermal fluctuations, in terms of a bath generated by the remaining vibrational modes which are not directly coupled to the charge-transfer process. This deserves a separate study. , and E z (t) = E Resonance + 0.05 (green). E Resonance is the electric field in which the greatest charge transfer between donor and acceptor is produced.