Quantum transport through open systems

The question of finding non-equilibrium stationary states, on which energy and particle transport are established, has been an important open question in Physics since very early stage of Physics. For classical systems, such as heat conductance and electric circuits, one can use phenomenological laws such as Fourier's law and Ohm's law, although even those laws are not derived directly from classical equation of motion. During the past many decades there have been many efforts on trying to deriving those laws from the first principle. When it comes to quantum systems, we do not even have such phenomenological laws. So the ultimate goal of this project is to develop such a theoretical framework to find non-equilibrium stationary states from the first-principle equation of motion of quantum systems.

In equilibrium statistical physics, an equilibrium state of a system (both classical and quantum) is given explicitly by the Boltzmann distribution, which requires only the Hamiltonian (and of course, its' eigenvalues and eigenvectors), bath temperature and chemical potential. For non-equilibrium states, even when we are given all those things, we still don't know how to write down explicitly the distribution function.

Therefore, what we are looking for is a blackbox accepting Hamiltonian, bath temperatures, chemical potentials, and if necessary coupling mechanism and coupling strength to baths, and giving a distribution function. In this project, we consider only quantum systems.

There are already some approaches to this problem, including: the Kubo formula, which treats the central systems as closed systems, the Landauer-Buttiker formula for non-interacting systems and the non-equilibrium Green's function method, which is capable of dealing with interaction. We decided to use another approach based on the open-system scenario: starting from the Schroedinger equation of the composite system of central system+bath, deriving an effective equation of motion of the central system only and then solving the effective equation. The baths are usually assumed to be in thermal equilibrium states all the time (i.e. baths have infinite degree of freedom and coupling strength is small so that back action from the central system to the bath can be neglected). Derivations of the effective equations usually involves additional approximation such as second-order perturbation and Markovian approximation. But those approximations can be relaxed and then more complicated effective equation can be derived. At the current stage of this project, we will not focus on deriving such equations, but rather on solving those equations, developing and applying this framework. One of such equations is called the Redfield equation. It makes the above two approximations but still makes use of the full dynamics of the central systems.

Note the whole picture mentioned above has already been suggested before we started this project, especially in discussion of relaxation towards equilibrium states. However, it is not widely used in describing non-equilibrium stationary states. We believe that it is because of the difficulty in solving the resulting open-system master equation such as the Redfield equation. Using direct methods, one needs to solve an eigenvalue problem of size 4N where N is the size of the system measured in qubits. If efficient methods of solving this equation can be found, then more applications on various physical models of this framework will for sure emerge.

The following lists several projects I have worked on.

BBGKY-like method for solving the open-system master equation

As we mentioned above, solving directly the density matrix involves linear system with size 4N. One way to reduce the dimension is to convert the equation of density matrices into equations of Green's functions. All single-particle physical quantities can be calculated from single-particle Green's function. So if we are interested in only such physical quantities, we only need to deal with N2 unknowns, {G(i,j), i,j=1...N}.

However, for interacting systems, unfortunately, those equations of single-particle Green's functions, denoted as G1, are coupled to two-particle Green's function G2. Furthermore, at the next order, equations of G2 are coupled to G3 and so on. This hierarchical structure exists too for equations of equilibrium Green's function. Solving the whole hierarchy is as hard as solve the original equation of density matrices. There has to be some ways to truncate and solve the hierarchy consistently and systematically. In equilibrium Green's function, the BBGKY (Bogoliubov-Born-Green-Kirkwood-Yvon) hierarchy and the cluster expansion provide exactly such a method.

Therefore, here we just need to develop a similar equation hierarchy for the Redfield equation. The only difference is at those additional bath operators in the Redfield equation: they are operators involving many-particle operators. One has to derive too a hierarchical form of those operators together with the interaction term in the Hamiltonian.

In its first-order form, we truncated at terms involving G2, neglected higher Gs and assumed

which is the expansion of G2 in terms of G1 according to Wick's theorem. A non-equilibrium Wick's theorem can be proved for non-interacting systems. Using this, the equation of G1 is closed. We solved this equation for G1 and compare it against exact solution from direct methods. Roughly speaking there is a order of a magnitude improvement of the accuracy from the zeroth order, which interaction is neglected totally, to the first-order, see in Figure Cluster1.

difference between G<sub>1</sub>s calculated from the BBGKY and the direct exact method difference between particle currents calculated from the BBGKY and the direct exact method

Figure Cluster1: (1) G1 calculated from the first-order BBGKY-like method is compared with those G1 calculated directly from the density matrices. dC,(0) stands for the comparison between G1 calculated from the zeroth-order (i.e. neglecting interaction) BBGKY-like method and those calculated directly from the density matrices. (2) Difference between particle currents calculated from two methods. Both figures demonstrate an essential improvement of accuracy.

In its second-order form, we truncated at terms involving G3, neglected higher Gs and assumed

G3p(-1)pG1G1G1 + Σp(-1)pG1G2
which is beyond Hartree-Fock approximation and is often used in cluster expansion. Using this, the equation of G1 and G2 is closed. Again we solved this equation and compare results against exact solution. The accuracy is improved further but more importantly, we now can find out non-trivial correlation in G2, see in Figure Cluster2.

difference between G<sub>1</sub>s calculated from the BBGKY and the direct exact method difference between G<sub>2</sub>s calculated from the BBGKY and the direct exact method difference between particle currents calculated from the BBGKY and the direct exact method

Figure Cluster2: (1) G1 calculated from the second-order BBGKY-like method is compared with those calculated directly from the density matrices. (2) Similar comparison on G2. (3) Difference between particle currents calculated from two methods. Clearly there are further improvement of accuracy from the first-order to the second-order BBGKY method. This is especially demonstrated in the figure of particle currents. In order to enlarge the difference between the first-order method and the exact one, large coupling strength has been used in those figure than the previous two.

Coherent-state representation approach

Using coherent-state representation, density matrices of N-qubit systems can be mapped to distribution functions over 2N complex variables and the equations of motion of the density matrices become differential equations of the distribution functions. This is usually done for evolution of pure dynamical system or system evolves towards equilibrium states.

In this work, we use the coherent-state representation approach to convert the 4N-dimensional Redfield equation into a stochastic differential equation with 2N complex variables, so called a generalized Fokker-Planck equation. Analytical expression of the non-equilibrium stationary states are derived for some systems. The accuracy of this method is around 6%, see in Figure Coherent. Manuscript in preparation.

For non-interacting systems, Green's function solved analytically and numerically in the coherent-state representation compared with those from the BBGKY-like method For interacting systems, Green's function solved numerically in the coherent-state representation compared with those from the BBGKY-like method

Figure Coherent: G1 calculated in the coherent-state representation is compared against those G1 calculated from the first-order BBGKY-like method for non-interacting bosonic systems in (1) and interacting systems in (2). For bosonic interacting systems, direct exact solution is not possible because the Hilbert spaces are infinite dimensional. dRe,NuB stands for the comparison between numerical solution in the coherent-state representation against the BBGKY-like method. The two figures show that the approach based on the coherent-state representation is working.

Kubo formula for open systems

We found in study of the Kubo formula for open systems that in order to study transport one has to take into account the coupling from the central system to the baths explicitly. In using the usual Kubo formula in transport studies, one assume the central system is a closed system and non-equilibrium properties are calculated from response of the equilibrium state to proper external variations. We found that first, the usual Kubo formula gives corrections towards the equilibrium state of the perturbed system not the presumed and widely believed non-equilibrium state; second, the usual Kubo formula is possible to be valid only when dealing with infinite systems, where boundary terms can always be neglected; and third, unphysical divergence always emerge from the "diagonal part" (or more generally, elements of density matrices between those states which have the same energy) in using the usual Kubo formula. Starting from the open-system Redfield equation instead of the close-system Schroedinger's equation solves all above problems. Results are found to be consistent with direct exact solutions.

Comparison among the first open Kubo formula, the second open Kubo formula and the exact solution

Figure open Kubo: D1 (D2) is the difference between the first (second) open Kubo formula and the exact solution. Nex, N1 and N2 are respectively particle number calculated from the exact solution, first and second open Kubo formula. Similarly particle currents Jex, J1 and J2 are defined. We see that the first open Kubo formula gives accurate results for all range of δ and on all physical quantities (D,N,J) while the second only generate good results on J. Density matrices from the second formula in fact diverge at all "diagonal" elements. But only for those physical quantities, which has zero "diagonal" elements, the second formula leads to proper values.

Thermal transport of spin chains

Using direct methods we studied thermal transport of spin chains and analyzed systems up to N=10. Connections between integrability and anomalous transport, which is widely believed by physicists and has been demonstrated by studies based on the usual Kubo formula, is challenged by our results. Instead we find a qualitative difference between non-interacting systems and interacting systems, which respectively in the case of spin systems correspond to Jz is zero and non-zero. This can be seen from the Jordan-Wigner transformation between spin and fermions. The same observation holds for tight binding chains with nearest neighbor interaction.

Temperature profile on spin chains spin chains with small Jz values Temperature profile on fermionic V-t chains

Figure Temperature: (a) Temperature profiles on spin chains. Notice there is a qualitative difference between zero and nonzero Jz. (b) Some extreme cases: very small Jz and small Bz. When Bz is too small, the idea of local energy and local equilibrium breaks down so it is hard to talk about temperature in that case. (c) Similar difference is observed on fermionic chains with nearest neighbor interaction.

Evolution of open systems from only the first principle?

The method based on the coherent-state representation is capable of dealing with very large systems and for their pure dynamical evolution, evolving towards thermal equilibrium and also towards non-equilibrium stationary states. Besides solving the open-system Redfield equation, there is much more we can do with this method. One of them is to investigate under what circumstances the composite system+bath configuration does drive the central system into stationary states, either equilibrium or non-equilibrium ones.

It is well known, and as we used and showed above, that if the bath is ideal that it has infinite degrees of freedom, it is initially in a thermal equilibrium state, and the coupling strength between the bath and the central system is infinitesimal, then the central system will be driven into a stationary state. In its not-so-rigorous derivation, usually second-order perturbation and Markovian approximation are used. So the question is how essential are all the above conditions and approximations for the central system to converge to the stationary states?

A well-known exactly solvable model of decoherence A sketch of a finite-size central system coupled to finite-size baths

Figure 1st_principle: (a) A well-known exactly solvable model of decoherence. (b) A central system is coupled to two finite-size baths. Approximations at different levels may be added into the description one by one. It is possible to use the coherent-state method to keep track of a pure dynamical evolution. In a sense this can be regarded as generalization of (a).

This might be answered if we can implement the pure dynamical evolution of the system+bath and add in various approximations and conditions one by one. For example, as in Fig 1st_principle we may take both the bath and the central system to be finite, of course with the bath being much large, and study their pure dynamical evolution with different initial condition of the bath's state. We can then set the initial bath states follows the Boltzmann distribution and check if the central system evolves towards stationary states when the bath size goes to larger and larger. From there we may also be able to understand the difference between Markovian and non-Markovian evolution.

In fact, answers to this question have deeper meaning than what they appears here. These answers could help us establish a framework of studies of general open system, including under what condition infinite-dimension assumption is valid or Markovian approximation is reasonable, or it is fine to assume the bath in equilibrium and so on. This question, how to study open systems, is really essential in quantum computing.

Comparison between this open-system approach and the NEGF method

In principle, this open-system approach and the NEGF method are trying to study the same physical question: transport on a composite system of system+bath. However, they deal with the same situation very differently.

It is so different that they even have different equilibrium limits. Let us first consider the NEGF configuration: lead+system+lead with equilibrium setup where both leads, which can be regarded as semi-infinite non-interacting chains, are in equilibrium states with the same temperatures and chemical potentials. One special case is when the central system is also non-interacting, say a finite non-interacting chain, which has the same natural with the leads. If we ask what is the stationary state of the whole system, one can easily imagine that it will the equilibrium state of the whole composite system of lead+system+lead. In fact, using Dyson's equation for Green's function, one can find exact solution for this case and confirm the above intuitive picture.

Let us now consider the open-system scenario: a central system is coupled to two baths, which may in fact as well be the same semi-infinite non-interacting chains but we start directly from the reduced equation of motion of the central system. From there, if the central system is a non-interacting one we find that the stationary state is the equilibrium state of the central system only.

What is the difference? In the case of NEGF, the equilibrium state is of the whole composite system while in the case of open-system scenario, the equilibrium state is of the central system only. Therefore, although they deal with the same physical question, there is no guarantee that they will have the same quantitative or even qualitative results. Then in what circumstance, will the two be directly comparable? Thinking about a two-spin Ising system for instance with coupling constant Jz. When there is no local magnetic fields on the spins, starting from the equilibrium state of the whole two-spin system,


where p(T) and q(T) depend on T but p(T)+q(T)=1/2, the reduced density matrix of one spin becomes

which has no T dependence at all. I learned from this example that in order to have the reduced density matrix to be approximately the induced equilibrium distribution from the composite system, local Hamiltonian such as the local magnetic field term should be an important part of the Hamiltonian of the central system. The second condition is quite intuitive that the coupling strength between the central system and bath should be very small compared to the Hamiltonian of the leads and the central system. In terms of the above example, the two conditions require

Bz >> Jz


Jz, intra >> Jz, inter

where Jz, inter refers to the coupling strength between the central system and the leads. Under those conditions, ideally the two methods should be comparable. In all of our previous example, coupling to baths are taken to be in the momentum space. We have derived some formulae for coupling in real space, using a configuration similar to Fig 1st_principle(b). Detailed analytical and numerical work has not yet started.