Including arbitrary geometric correlations into one-dimensional time-dependent Schr\"{o}dinger equations

The so-called Born-Huang ansatz is a fundamental tool in the context of ab-initio molecular dynamics, viz., it allows to effectively separate fast and slow degrees of freedom and thus treating electrons and nuclei at different mathematical footings. Here we consider the use of a Born-Huang-like expansion of the three-dimensional time-dependent Schr\"odinger equation to separate transport and confinement degrees of freedom in electron transport problems that involve geometrical constrictions. The resulting scheme consists of an eigenstate problem for the confinement degrees of freedom (in the transverse direction) whose solution constitutes the input for the propagation of a set of coupled one-dimensional equations of motion for the transport degree of freedom (in the longitudinal direction). This technique achieves quantitative accuracy using an order less computational resources than the full dimensional simulation for a prototypical two-dimensional constriction.


Introduction
Nanoscale constrictions (sometimes referred to as point contacts or nanojunctions) are unique objects for the generation and investigation of ballistic electron transport in solids. Studies of such systems have been inspired by the pioneering investigations by Sharvin in the mid-1960s [1]. Today, advances in the fabrication techniques like direct growth of branched nanostructures [2], electron beam irradiation [3], thermal and electrical welding [4] or atomic force microscope [5] have allowed to control the size and composition of nanojunctions for creating devices with desired functionalities. In this respect, a number of nanodevices based on nonjunctions like the single electron transistors [6,7], field effect transistors [8,9] and heterostructure nanowires [10,11] have been recently reported which promise great performance in terms of miniaturization and power consumption.
In the design of these nanostructures, simulation tools constitute a valuable alternative to the expensive and time-consuming test-and-error experimental procedure. Today, a number of quantum electron transport simulators are available to the scientific community [12][13][14][15][16]. The amount of information that these simulators can provide, however, is mainly restricted to the stationary regime and therefore their predicting capabilities are still far from those of the traditional Monte Carlo solution of the semi-classical Boltzmann transport equation [17]. This limitation poses a serious problem in the near future as electron

Single-electron time-dependent Schrödinger equation in a Born-Huang-like basis expansion
As we have explained in the introduction, it is our goal to reduce the computational burden associated to the solution of an ensemble of effective single-electron 3D SSE [32]. Therefore, we consider our starting point to be the 3D TDSE of a single (spin-less) electron in the position basis, i.e.: i ∂ ∂t Ψ(x, y, z, t) = H(x, y, z)Ψ(x, y, z, t), where we have used atomic units, and x, y and z represent the three spatial coordinates. In Equation 1, H(x, y, z) is the full Hamiltonian of the system: H(x, y, z) = T x + T y + T z + V(x) + W(x, y, z), which has been assumed to be time-independent for simplicity. The time-dependence on the scalar potentials V(x) and W(x, y, z) will be discussed in later sections. In Equation 2, T x = − 1 2 ∂ 2 ∂x 2 and V(x) are, respectively, the kinetic energy and the scalar potential associated to the longitudinal degree of freedom x, while T y = − 1 2 ∂ 2 ∂y 2 and T z = − 1 2 ∂ 2 ∂z 2 are the kinetic energies associated to the transversal degrees of freedom y and z. The scalar potential W(x, y, z) includes any other scalar potential that is not purely longitudinal, which is responsible of making the solution of Equation 1 non-separable.
It is convenient at this point to rewrite the Hamiltonian in Equation 2 in terms of longitudinal and transverse components as: where H ⊥ x (y, z) is the transverse Hamiltonian defined as: An eigenvalue equation associated to the transverse Hamiltonian can now be introduced as follows: where E k (x) and φ k x (y, z) are the corresponding eigenvalues and eigenstates respectively. The eigenstates φ k x (y, z) form a complete basis in which to expand the Hilbert space spanned by the variables x, y, and z. Therefore, the 3D wavefunction in Equation 1 can be expressed in terms of transverse eigenstates φ k x (y, z) as: where χ k (x, t) = dydzφ k x (y, z)Ψ(x, y, z, t) are complex longitudinal coefficients associated to the transverse eigenstate φ k x (y, z). Unless otherwise stated all integrals are evaluated from −∞ to ∞. It is important to note that since the longitudinal variable x appears as a parameter in Equation 5, the transverse eigenstates obey the following partial normalization condition: In addition, the longitudinal complex coefficients χ k (x, t) in Equation 6 fulfill, by construction, the condition: The wavefunction expansion in Equation 6 can now be introduced into Equation 1 to obtain an equation of motion for the coefficients χ k (x, t) (see Appendix A): where E k (x) are effective potential-energies (that correspond to the eigenvalues in Equation 5) and F kl (x) and S kl (x) are geometric (first and second order) coupling terms, which read: Since the transverse eigenstates φ k x (y, z) are real, the term F kk is zero by construction. The other terms in Equation 10 dictate the transfer of probability presence between different longitudinal coefficients χ k (x, t) and, therefore, will be called geometric non-adiabatic couplings (GNACs). Accordingly, one can distinguish between two different dynamics regimes in Equation 9: (i) Geometric adiabatic regime: it is the regime where F kl and S kl are both negligible. Thus, the solution of Equation 9 can be greatly simplified because it involves only one transverse eigenstate. (ii) Geometric non-adiabatic regime: it is the regime where either or both F kl and S kl are important.
Thus, the solution of Equation 9 involves the coupling between different longitudinal coefficients and hence more than one transverse eigenstate.
Interestingly, the prevalence of the regimes (i) or (ii) can be estimated by rewriting the first order coupling terms F kl (x) as (see Appendix B for an explicit derivation): That is, the importance of non-adiabatic transitions between transverse eigenstates depends on the interplay between the transverse potential-energy differences E l (x) − E k (x) and the magnitude of the classical force field ∝ ∂ ∂x W(x, y, z). The geometric adiabatic regime (i) is reached either when the classical force field is very small or the energy differences E l (x) − E k (x) are large enough. In the adiabatic regime only the diagonal terms, S kk , are retained, which induce a global shift of the potential-energies E k (x) felt by the longitudinal coefficients χ k (x, t). In this approximation, the longitudinal degree of freedom moves in the potential-energy provided by a single transverse state, E k (x). This regime is analogous to the so-called Born-Oppenheimer approximation in the context of molecular dynamics [43], where the term S kk is often called Born-Oppenheimer diagonal correction [44]. As it will be shown in our numerical example, the evolution of the system can be governed either by the geometric adiabatic or nonadiabatic regime depending on the particular spatial region where the dynamics is occurring.
Let us notice at this point that the time-dependence of the Hamiltonian in Equation 2 may come either due to a purely longitudinal time-dependent scalar potential V(x, t) or through the time-dependence of the non-separable potential W(x, y, z, t). If the time dependence is added only through V(x, t), then nothing changes in the above development. Contrarily, if a time-dependence is included in W(x, y, z, t), then the eigenstate problem in Equation 5 changes with time and so do the effective potential-energies E k (x, t) and the first and second order GNACs F kl (x, t) and S kl (x, t). As it will be shown later, in this circumstance, Equation 5 should be solved self-consistently with Equation 9.
Before we move to a practical example implementing the above formulation of the 3D TDSE, let us emphasize that it is the main goal of the set of coupled equations in Equation 9 to allow the evaluation of relevant observables in terms of 1D wavefunctions only. In this respect, let us take, for example, the case of the reduced probability density ρ(x, t) = dydzΨ * (x, y, z, t)Ψ(x, y, z, t). Using the basis expansion in Equation 6, ρ(x, t) can be written as: and using the condition in Equation 7 the above expression reduces to: Therefore, according to Equation 13, the reduced (longitudinal) density is simply the sum of the absolute squared value of the longitudinal coefficients Similarly, other relevant observables, such as the energy can easily be derived using the expansion in Equation 6 (see Appendix C).

Application of the method to a prototypical constriction
The above formulation can be cast in the form of a numerical scheme to solve the 3D TDSE, which we will call, hereafter, geometrically correlated 1D TDSE (or in brief GC-TDSE). The scheme can be divided into two different parts corresponding to their distinct mathematical nature. The first part involves the solution of the eigenvalue problem in Equation 5 which allows to evaluate the transverse eigenstates φ k x (y, z) and eigenvalues E k (x) as well as geometric non-adiabatic couplings F kl (x) and S kl (x) in Equation 10. These quantities are required in the second part of the algorithm to solve the equation of motion of the longitudinal coefficients χ k (x, t) in Equation 9, which ultimately allow us to evaluate the observables of interest.
In what follows we discuss these two aspects of algorithm for a prototypical 2D geometric constriction whose geometry does not change in time. We consider one degree of freedom in the transport direction and one degree of freedom in the transverse (or confinement) direction. The generalization to a 3D system, i.e., with two transverse degrees of freedom, is straightforward and does not add any physical insight to the 2D case. As it will be shown, the transverse eigenstates and eigenvalues as well as the geometric non-adiabatic couplings are, for a time-independent constriction, functions that can be computed once and for all. That is, the effective potential-energies E k (x) and the non-adiabatic couplings F kl and S kl are computed only at the very beginning of the GC-TDSE propagation scheme. For more general time-dependent constrictions, possibly with no analytical form of W(x, y, z) and hence of φ k (y) and E k (x), the only change in the algorithm is that Equation 5 has to be solved, self-consistently, together with Equation 9 at each time step.

Evaluation of transverse eigen-states (and values) and geometric non-adiabatic couplings
Let us consider the case of a 2D nanojunction represented by the scalar potential: where L 1 (x) and L 2 (x) define the shape of the constriction. Given the 2D Hamiltonian, where H ⊥ x (y) = T y + W(x, y). The wavefunction for a 2D constriction in terms of Born-Huang expansion can be written as, , and the transverse states φ k x (y) are solutions of a free particle in a 1D box whose length depends upon the longitudinal variable x, i.e.: The associated eigenvalues are given by: where we have defined . These energies, parametrically dependent on the longitudinal variable x, define the effective potential-energies where the coefficients χ k (x, t) evolve on.
To evaluate the first and second order coupling terms F kl (x) and S kl (x), we need to rely on a particular form of the constriction. Depending on the specific form of L 1 (x) and L 2 (x), different geometrically bounded constrictions can be conceived (see for example the panels (a) and (c) of Figure 1). Given the states in Equation 16 and a particular shape of the constriction (defined in Equation A17 of Appendix D), it is then easy to evaluate the non-adiabatic coupling terms F kl and S kl (see panels (b) and (d) of Figure 1).
The two different constrictions in Figure 1 serve well to gain some insight into the general form and dependence of the effective potential-energies E k (x) as well as of the GNACs in Equation 10. Geometries changing more abruptly lead to sharper effective potential-energies E k (x) and more peaked non-adiabatic coupling terms F kl (x) and S kl (x). Sharper constrictions are thus expected to cause larger non-adiabatic transitions and hence to involve a larger number of transverse eigenstates requiring a larger number of longitudinal coefficients in order to reconstruct the reduced (longitudinal) density in Equation 12. On the contrary, smoother constrictions should yield softer non-adiabatic transitions and hence involve a smaller number of transverse eigenstates and equivalently a smaller number of longitudinal coefficients.

Time-dependent propagation of the longitudinal coefficients
Given the effective potential-energies E k (x) and the non-adiabatic couplings F kl (x) and S kl (x), one can then easily find a solution for the longitudinal coefficients in Equation 9. Here, we consider the dynamics of an electron that impinges upon a constriction defined by Equations 14 and A17 in Appendix D using the particular set of parameters, {A, B, a 1 , a 2 , γ} = {180, 220, 630, 870, 20}, which corresponds to panel (c) of Figure 1. For this particular set of parameters, the symmetry of the states defined in Equation 16 forbids transitions between odd and even states, i.e., , the coupling between odd and even states is zero, i.e., F kl (x) = S kl (x) = 0, ∀k + l = odd.
We will consider two different initial states. On one hand, the initial wavefunction Ψ(x, y, 0) will be described by: where φ 1 x (y) is the transverse ground state defined in Equation 16, and ψ( x is a minimum uncertainty (Gaussian) wavepacket with initial momentum and dispersion k 0 = 0.086 a.u. and σ x = 80 a.u. respectively, and centered at x 0 = 300 a.u. (while N is a normalization constant). On the other hand, we will consider the initial state to be defined by: where now both ξ(y) and ψ(x) (defined above) are Gaussian wavepackets. In particular, ξ(y) = M exp −(y−y 0 ) 2 2σ 2 with y 0 = 200 a.u., σ y = 20 a.u. (and M a normalization constant). The probability densities |Ψ(x, y, 0)| 2 associated to the above two initial states can be seen, respectively, in panels (a) and (b) of Fig. 2. Given the initial states in Equation 19 and Equation 20, we can then deduce the corresponding longitudinal coefficients as follows: While the initial state in Equation 19 corresponds to χ k (x, 0) = δ k1 ψ(x), the state in Equation 20 involves a larger number of transverse eigenstates. Note that this second initial condition may be more realistic in In addition, we also show the absolute squared value of the longitudinal coefficients, i.e., |χ k (x, t)| 2 , evaluated using the GC-TDSE. Alternatively, in the right panels of Figs. 3 and 4 we plot the population of each transverse state, as a function of time using the GC-TDSE. The initial state in Equation 19 yields P k (0) = δ k1 . This can be seen in the right hand panel of Fig. 3. This value stays constant until the wavepacket hits the constriction at around t = 2500 a.u. At this moment, non-adiabatic transitions between different transverse states start to occur and lead to complicated interference patterns at later times (see, e.g., the reduced density ρ(x, t) at t = 5010 a.u.). The number of significantly populated transverse states increases up to six (while up to eleven states are required to reproduce the exact reduced density up to a 0.1% error). Among these states, only odd transverse states are accessible due to the symmetry of the initial state (as we noted in Equation 18). Since the mean energy of the initial state in Equation 19 ( Ê = 0.0037 a.u.) is higher than the barrier height of the first effective potential-energy in Figure 1 (d) (max(E 1 x ) = 0.0028 a.u.), one could naively expect a complete transmission of the wavepacket χ 1 (x, t). However, due to the effect of the non-adiabatic coupling terms F kl (x) and S kl (x), χ 1 (x, t) looses a major part of its population in favour of higher energy transverse components that are reflected by much higher effective potential-energy barriers.
Starting with the second initial state in Equation 20 there are up to seven transverse states populated at the initial time, all of which are again odd states due to the symmetry of the initial conditions (see Figure  4(d)). Once the wavepacket hits the constriction at t = 3000 a.u, states 3, 7 and 9 become more dominant than the previously dominant states 1, 3, and 5. Overall, up to 15 states become important to reproduce the exact reduced longitudinal density within a 0.1% error. Given the above two examples, it seems clear that the specific form of the impinging wavefunction does not play a determinant role in the scaling of the number of relevant transverse states even though it do effect the total number of states N e required to numerically evaluate Equation 9.
where x m is the center of the nanojunction in the longitudinal direction (i.e., with numerical value 750 a.u.) and t f is the time at which no probability density (i.e., less than 0.1%) remains inside the constriction. In Figure 5 we show results for applied bias 0.0075 a.u. ≤ V ext ≤ 0.045 a.u. and two different number of states N e = 15 and N e = 25. As expected, higher applied bias lead to a more vigorous collision of the  wavepacket against the constriction due to a higher longitudinal momentum/energy, which allows higher energy transverse states to be populated.

General discussion
The GC-TDSE algorithm discussed in the previous sections has a clear computational advantage over the solution of the full 2D TDSE. This is particularly so when the quantities E k (x), F kl (x) and S kl (x), involved in the equation of motion of the longitudinal coefficients χ k (x, t), are time-independent functions. For a time-independent transverse Hamiltonian, the quantities E k (x), F kl (x) and S kl (x) are computed only once before the propagation of the longitudinal coefficients and thus the computational cost of the GC-TDSE resides, mainly, on the propagation of the 1D longitudinal coefficients.
Let us provide some numbers to get an estimate of the numerical efficiency of the GC-TDSE algorithm. Consider the numerical solution of the full 2D TDSE in a grid. For a number of grid points {n x , n y } = {1500, 400}, the resulting Hamiltonian has a dimension (n x × n y ) 2 . Alternatively, the size of the Hamiltonian involved in the propagation of the longitudinal coefficients of the GC-TDSE algorithm is (n x × N e ) 2 , where N e is the number of transverse eigenstates. One can then estimate the numerical efficiency of one method over the other by simply evaluating the ratio (n x × n y ) 2 /(n x × N e ) 2 . Thus, for time-independent transverse potentials W(x, y, z), the computational reduction associated to the GC-TDSE is n 2 y /N 2 e . Note that the benefits of the GC-TDSE would be even more noticeable when applied to a 3D problem, for which the above ratio would become (n 2 y × n 2 z )/N 2 e . As we have seen in the above section, the number of required transverse states N e is a function of the abruptness/smoothness of the constriction, but also of the energy of the impinging wavepacket. Therefore the computational advantage of the GC-TDSE method over the full dimensional TDSE is clearly system-dependent. Slow wavepackets impinging upon smooth constrictions would maximize the benefits of the GC-TDSE. Contrarily, very energetic electrons colliding against abrupt constrictions would certainly minimize its benefits. In this respect, we must note that the GNACs have a clear dependence upon the profile of the constriction. In particular, the second order coupling terms S kl (x) will be sharply peaked for very abrupt constrictions (see the important differences in the size and sharpness of the S kl (x) in Figure 1 for two different constrictions). Therefore, due to the non-unitary character of the equations of motion of the longitudinal coefficients, very abrupt constrictions may demand very fine grids in practice.
Finally, let us mention that whenever the transverse Hamiltonian in Equation 4 is time-dependent, the advantage of the GC-TDSE method compared to the solution of the full dimensional TDSE is not so obvious. As we have already noticed, for a time-dependent transverse potential W(x, y, z, t), the eigenvalue problem in Equation 5 must be solved self-consistently with Equation 9, i.e, at each time step. Then, a comparison of the GC-TDSE and the full dimensional TDSE in terms of numerical efficiency will depend on the specific performance of the eigensolver utilized to evaluate the transverse eigenvalues, E k (x, t), and eigenstates φ k x (y, z).

Conclusions
In this work we have proposed a new method, named GC-TDSE, that allows to include arbitrary geometric correlations between traversal and longitudinal degrees of freedom into a coupled set of 1D TDSE. Our motivation for the development of this method was, initially, a further reduction of the dimensionality of the 3D Schrödinger-like equations that result from a (Monte Carlo) SSE approach to quantum electron transport in open systems (valid for Markovian and non-Markovian regimes) that we have recently proposed [32]. Nevertheless, the method presented here is general and allows to reduce the dimensionality of quantum systems with geometrical correlations among different degrees of freedom, which could be of great utility in different research fields.
For smooth time-independent constriction profiles under low applied bias, our GC-TDSE method implies up to three orders of magnitude less computational resources than solving the full 3D TDSE directly. For very high applied bias or time-dependent constrictions profiles, the GC-TDSE may still be significantly cheaper than the solution of the full 3D TDSE, but would require introducing approximations to the solution of the potential-energies E k (x, t) and the GNACs (F kl (x, t) and S kl (x, t)). We thus expect the GC-TDSE presented here to trigger future investigation for making it robust against more extreme electron transport conditions and to inspire new ways of looking at the many-body problem.

Conflicts of Interest:
The authors declare no conflict of interest. φ l |H ⊥ x |φ k = E k φ l |φ k + φ l |(H ⊥ ) |φ k + E l φ l |φ k (A9) Using the relation in Equation A7 in Equation A9 we get, In the above equation we have made use of the fact that φ l |H ⊥ |φ k = (E k φ l |φ k ) = 0, when k = l Therefore, Which can be written in the position representation and using the same nomenclature as used in the main text as follows, Using Equation 4 it is easy to see that only the potential function W(x, y, z) will survive after the partial derivative with respect to the variable x. so we can equivalently write Equation A12 as,