Discrete Numerical Model

To explore the eigenvalue spectrum of the present open-system model and those which will be investigated later, let us consider a finite-difference approximation to the Liouville equation (2.3) which reduces to a finite matrix whose eigenvalues may be readily computed. Let me emphasize that only the spatial coordinates will be discretized; time remains continuous so that the partial differential (and eventually integro-differential) Liouville equation will be reduced to a set of coupled ordinary differential equations with respect to time.

This particular situation requires some discussion. Throughout the computational physics literature, discussions of stability always involve a discretization with respect to time. Because one is accustomed to dealing with continuum equations whose behavior is known to be stable (or at least physical) the common assumption that any instability must be a result of the discretization scheme is generally correct. However, a different situation is being studied here. The validity of the equations themselves (or more precisely the boundary conditions) is the issue. If a discrete-space, continuous-time model is unstable, there will be no time discretization which will correct this instability. On the other hand, we wish to assume that the stability of the discrete-space, continuous time model will be indicative of the stability of a continuous-space, continuous-time model. As mentioned before, this connection requires a logical induction.

The position coordinates x will be taken to be elements of a uniformly spaced mesh: . The dependent quantities such as the wavefunction and density matrix then take on discrete values also, which will be denoted by , and . Using the simple finite-difference approximation , the Hamiltonian (2.2) becomes

for i, j not on one of the boundaries. To incorporate the boundary conditions, it is best to think of adding an additional mesh point at each end of the domain (points and ), and specifying the value of the wavefunction on those points. For example, to apply the homogeneous Dirichlet conditions for a particle in a box, we would set and . Inserting these conditions into (3.24) completely defines the matrix for . Similarly, if we wanted to apply Neumann conditions , we would set .

Writing the Liouville equation (2.3) on the finite-difference basis gives

where the tetradic nature of is made explicit. This discrete representation of may be derived from (3.24) and is

Again, the elements adjacent to a boundary require special attention.

To evaluate the eigenvalues of and other superoperators, we must map the tetradic onto an ordinary matrix, so that conventional eigenvalue algorithms may be applied. To do so for the finite, discrete case, we may map the density matrix onto a singly subscripted vector of dimension by with . Note that with this mapping the inner product between two operators (2.4) becomes the ordinary inner product between two vectors. The mapping of the tetradic onto an matrix follows immediately. The matrix representing was actually constructed for N = 8 (resulting in a matrix for ) using the potential illustrated in Fig. 4.

Figure 4. Potential used in evaluating eigenvalue spectra of Liouville superoperators in the discrete model. This potential was chosen to have both a driving field and a barrier.

Figure 5. Eigenvalue spectra of the Liouville operator for a small model closed system with the potential shown in Fig. 4. If the system is taken to be conservative, the resulting eigenvalue spectrum is shown in (a). All eigenvalues are purely real, as expected. In (b) a damping term has been added, leading to negative imaginary parts for most eigenvalues.

In the second case the model system is taken to be closed, but damped. The Fokker-Planck damping operator (3.23) may be written in discretized form as

This form preserves the important properties of . To illustrate the effect of dissipation on the spectrum of , the zero temperature limit () was taken (so that the first term, describing fluctuations, vanishes) and the damping constant (where ) was used. The resulting eigenvalue spectrum is shown in Fig. 5(b). Negative imaginary parts have been introduced into all the eigenvalues (except possibly one eigenvalue which is equal to zero within the numerical roundoff error, which presumably represents the ground state). These negative imaginary parts lead to damped motion, as expected.

Now with this background we can consider the case of the open-system boundary conditions (3.19) (zero diagonal gradient). The simplest finite-difference approximation for the condition (3.19) is

for i or j equal to 1 or N. Thus, the open-system Liouville superoperator (for open system, reversible) is obtained by inserting boundary values and (and the expressions obtained by transposing the indices) into (3.26). For the sake of completeness, let us write down the elements of which are affected by the boundary conditions:

The non-Hermiticity of follows from these expressions. For example, , but . The boundary conditions have caused elements of to be canceled in a way which breaks the Hermitian symmetry. The resulting eigenvalue spectrum is plotted in Fig. 6(a). The non-Hermiticity of leads to some eigenvalues with nonzero imaginary parts. It is apparent that these eigenvalues occur in complex conjugate pairs, with both positive and negative imaginary parts present. This is a consequence of the time-reversal symmetry of both the Liouville equation and the open-system boundary conditions (3.19). The eigenvalues with positive imaginary parts produce growing exponential solutions to the Liouville equation, which would prevent any approach to steady state. This open-system model is thus physically unacceptable.

Figure 6. Eigenvalue spectra for open systems using the boundary conditions of equation (3.19). If the boundary conditions are changed so as to open the system, nonzero imaginary parts are generated, as in (a). Because the boundary conditions are time-reversible, these imaginary parts occur in conjugate pairs. If a damping term is added as in (b), most, but not all, imaginary parts are negative. The few eigenvalues with positive imaginary parts are sufficient to render the model unstable. Stability can be achieved by increasing the damping rate, leading to the spectrum (c).

One might speculate that the problem of growing solutions could be due to the absence of damping in the model. To test this, let us add in the Fokker-Planck damping term (3.27), as we did for the closed-system model. With the same damping constant [] as before, the resulting eigenvalue spectrum for is that shown in Fig. 6(b). The addition of damping clearly does not solve the stability problem because it does not remove the positive imaginary parts. In fact, a larger damping constant does lead to a stable model, as shown in Fig. 6(c), where was used. All the eigenvalues now have negative imaginary parts, except for a doubly degenerate eigenvalue at zero (which must be present because of the invariance of and ).

Thus, modeling an open system by applying the boundary conditions (3.19) will work only if the rate of damping within the system is sufficiently large (or, for the case of electron transport, if the mobility is sufficiently low). The minimum acceptable damping rate depends upon the magnitude of the imaginary parts of the eigenvalues of for the undamped system, which in turn depends upon the form of the potential. In fact, the potential of Fig. 4 was chosen because it produces larger imaginary parts than potentials with greater symmetry. All this adds up to a very unsatisfactory formulation for an open-system model. The problems may be traced to the time-reversal symmetry of the boundary conditions. To obtain a proper formulation, this symmetry must be broken.