Introduction

Modern electronics relies on the ability of devices to switch internal states in a controllable manner. With the recent advantages in manufacturing complex nanoscale systems with engineered wavefunctions and energy surfaces, prediction of characteristic timescales has become an important and challenging task for a theory1,2. In particular, systems such as metal/oxide/metal memory cells3, spintronic devices4, including quantum dots of single atoms and molecules on metallic surfaces5,6,7 and organic solar cells8 are of a quantum nature and exhibit a complex energy landscape, resulting in a dynamics with multiple characteristic timescales.

In general, multiscale transition dynamics appear in any system with a complex energy surface and are enhanced when the system interacts with an external reservoir. These dynamics manifest themselves in the presence of metastable states, in which the system stays for a long time. The existence of multiscale transition dynamics is not limited to a particular class of compounds and is possible in a wide range of systems, including chemical reactions9, polarons in molecular junctions10,11 and black holes12. In a classical picture, metastability occurs for any isolated local energy minima. To describe quantum objects, one should consider metastable states as those from which a transition to lower energy states is prohibited by selection rules. The destruction of metastable states is possible via thermal activation to higher-lying energetic configurations or via tunneling.

Nanoscale systems are of a quantum nature and the transitions occur between different multi-electron configurations. A natural way to observe metastability is by studying the relaxation dynamics in the presence of a weak coupling to the environment. It is usually represented as a reservoir of non-interacting particles and is described by the Hamiltonian with states k forming a band. This reservoir is coupled to the system by particle exchange reflected in a hybridization Hamiltonian where subscript α enumerates localized states of the system and coefficients γαk represent their hybridization amplitudes. It is thus of crucial importance to provide a microscopic description of the dynamics which takes into account the full many-particle nature of the eigenstates of the systems, decoherence effects and charge/energy transfer due to the coupling with the surrounding environment.

Various involved computational approaches such as time-dependent numerical renormalization group13, time-dependent density matrix renormalization group14,15,16,17, time-dependent functional renormalization group18,19,20 and other methods based on quantum Monte Carlo algorithms21,22,23 do not allow one to access the behavior of the system at large times. On the other hand, the reduced density matrix theory24,25,26,27,28,29 provides a feasible way to access long-time properties of an interacting open quantum system. An exact version of this formalism can be coupled with the numerical methods mentioned above to obtain long-time dynamics at certain limits30,31. A frequently used approximation is described by the Redfield equation32,33 which has been used in many contexts including nuclear magnetic resonance34,35,36, quantum optics37,38,39 and tunneling through quantum dots40,41,42,43.

In the reduced density matrix theory, the environment acts as a reservoir of particles and energy for the system and is often assumed to be of Markovian nature. A Fermi golden rule approach44,45 allows one to write the rate equations for diagonal elements of the density matrix ρS(t) = trEρ(t). The method of full counting statistics46,47,48 builds upon it to allow for an efficient evaluation of time-dependent observables and correlation functions. This approach neglects quantum interference effects between different eigenstates. To take them into account, one can employ the Lindblad formalism49,50,51. For a range of non-interacting problems the standard choice of transition operators as local creation and annihilation operators is known to be valid52,53,54. The Lindblad formalism is also widely used for correlated systems interacting with an empty bath. This is a typical setup in the field of quantum optics where driven quantum systems emit photons to the vacuum but does not receive them back39. Much less is known on the applicability of the Lindblad approach for the description of the correlated quantum systems that can exchange particles with the bath. This situation is relevant for the description of nanoclusters, molecules, quantum dots, etc. interacting with a substrate.

In this report we consider the correlated quantum triangular cluster exchanging electrons with the fermionic reservoir. We construct an extension of the standard Lindblad one-particle formalism for correlated open quantum systems exchanging particles with the bath. Using the obtained formalism, we describe the physics of metastability in the relaxation process of the correlated system.

We follow the standard approach of constructing the master equation for the reduced density matrix by employing a perturbation theory in powers of HI and considering the terms up to the second order (so-called a sequential tunneling approximation40,42). The main distinction of the approach from the standard schemes that the transitions operators are constructed as matrices in the space of many-particle states. The time evolution in a leading contribution involves single-electron exchange with a reservoir. At low temperatures this results in a metastable many-body state of the system. We found that the metastablity arises from the specific structure of many-body spectrum and the coupling to a bath via the one-particle exchange channel. We consider a cluster of three atoms or quantum dots with local Coulomb interaction, which has a rich multiplet structure55,56 due to the interplay between frustration and interaction. The parameters of the cluster can be changed by using different isotopologues of nanoclusters or molecules57,58, while the substrate acts as a reservoir. The quantum simulation using ultracold atoms in optical traps is accessible. One can change the geometry of clusters and an interaction parameters directly59 and the environment is represented by a surrounding cloud of particles.

Results and Discussion

We considered the relaxation dynamics of an open quantum system. It is described by the master equation

The first Liouvillean term governs the unitary evolution of the system with the Hamiltonian . The second one takes into account the influence of the environment. It is obtained by applying a sequential tunneling approximation as described in Methods section.

We found that splits into so called one-particle and many-particle contributions. The first one is represented in terms of single-particle operators , cα as in the well-known Lindblad formalism approach39. The second one has an essentially many-particle nature and vanishes in the following cases: 1) infinite temperature T = ∞, 2) infinite chemical potential μ = −∞ or μ = +∞. This implies that the influence of correlations in the system on its relaxation dynamics comes into play only when both the temperature and the chemical potential of the reservoir are not infinite. The many-particle part of the superoperator is responsible for the relaxation to the Gibbs state at large time scales with the temperature T and chemical potential μ determined by the reservoir. This is also the case in a Fermi golden rule approach, but in contrast to it, master equation formalism describes the evolution of both diagonal and off-diagonal elements of density matrix ρS.

We have modeled the multiscale dynamics in an asymmetric triangular quantum dot with local Coulomb interaction which is shown in Figure 1. The expression for the Hamiltonian of the system under investigation is

Figure 1
figure 1

Schematic representation of the triangular cluster with one site coupled to the reservoir with temperature T and chemical potential μ.

It is modeled by the Fermi-Hubbard Hamiltonian in the half-filling regime with on-site interaction strength U and hopping amplitudes V = 1, V′ = 2. Coupling strength is characterized by a parameter J.

Full size image

Here U is the on-site interaction strength, V, V′ = 2V are hopping amplitudes. This setup represents a minimal frustrated interacting model, where the frustration occurs due to an antiferromagnetic exchange interaction in the chain of odd number of nodes. The interplay between frustration and interaction gives rise to a nontrivial structure of the energy spectrum which can be tuned by parameters of the model. This allows us to probe different relaxation regimes by changing the single parameter - interaction strength U. We choose hopping amplitude V as an energy scale of the system and break the C3 rotational symmetry by setting V′ ≠ V. This separates dynamical effects originating from electronic correlations from those specified by geometry. States of the system are classified using integrals of motion – the number of electrons N in the cluster, the square of its spin momentum S2 and the z-component Sz of its spin momentum. Note that the chemical potential on the cluster is μC = U/2. The triangular cluster lacks the electron-hole symmetry and its ground state does not necessarily correspond to N = 3 electrons. The low energy part of the spectrum of the Hamiltonian (2) is shown in Figure 2 as a function of U. At U ≈ 2.64 the ground state changes from N = 2 triplet state to N = 3 duplet state. There are also level crossings of excited states, which may result in distinctly different decay behavior when the system is brought to the contact with the bath. We consider dynamics when the bath is weakly coupled to the cluster. The details of the formalism used for the description of the relaxation processes of the open system in the regime of strong correlations are presented in the Supplementary Materials.

Figure 2
figure 2

Lowest energy terms as functions of the on-site interaction strength U.

Different colors are used for terms with different occupation numbers, line thickness indicates their multiplicity 2S + 1. Level crossings A1, A2 and A3 show points where tuning U leads to the change in dynamics behavior.

Full size image

The relaxation process is determined by the selection rules, ΔN = ±1, ΔSz = ±1/2. These are imposed on transitions between many-body states because the elementary process generated by the interaction Hamiltonian is a transfer of one electron. The resulting metastability includes (i) slow mixing of degenerate ground states which are not coupled by direct transitions and (ii) slow decay of non-ground states from which direct transitions to lower energy ones are not allowed. In both cases the relaxation happens through transitions to intermediate states lying higher in energy and its rate R1 ∝ Jexp(−ΔE/T) depends on the temperature of the reservoir and the effective coupling parameter J. The increase in T causes the system to escape faster from the metastable state. This holds in the time region where the thermal activation is the only mechanism of metastability decay. For larger times the effects of next order perturbations come into play. They allow exchanges with the reservoir by two particles at a rate R2 ∝ J2/V. Transitions between states with the same occupation number and ΔSz = ±1 act as another mechanism of the decay process and should be taken into account at low temperatures when the thermal activation rate is negligible.

We demonstrate discussed features of relaxation in the decay of a fully spin polarized state N = 3, S = 3/2, Sz = 3/2 at small temperature. It can be prepared by applying an external static magnetic field to the system. The behavior of dynamics depends on the interaction strength, so different values of U are considered. They are taken from the regions delimited by level crossing points A1, A2 and A3 as shown in Figure 2. The dynamics are manifested in the expectation values of observables Sz(t) and N(t) as shown in Figure 3. Their time dependence levels out at characteristic times Jt ~ 10–100. For U = 2.0 (lying below A1) there is a transition channel from the initial state to the ground state, so that metastability does not appear. For U = 2.64 (lying at A1) the system relaxes to an equiprobable mixture of a singlet state and two doublet states with 〈N〉 = 8/3 and 〈Sz〉 = 0 at large times. This mixing becomes metastable in the region above A1. The system moves to the state with Sz = 1/2 of the ground doublet while the state with Sz = −1/2 can only be reached through an excitation to the singlet state. This process might be accelerated by tuning the system closer to A1 – compare U = 3.0 and U = 3.5. Passing through A2 point changes the number of transition channels leading to the ground state from 2 to 1. The curve for U = 4.5 levels out more slowly than the curve for U = 3.5 because the intermediate state N = 2, S = 1 becomes higher in energy than the initial state. Increasing U above A3 point blocks the transition through the second intermediate state N = 4, S = 1 as well. For U = 6.0 an initial fully spin polarized state becomes metastable with respect to the decay.

Figure 3
figure 3

The decay of a fully spin polarized state in the frustrated triangular cluster coupled to the reservoir at temperature T = 0.02.

The interaction parameter U is tuned as indicated by different colors. Dynamics of (a) cluster magnetization 〈Sz〉 and (b) cluster occupation 〈N〉 is presented. One can see that tuning U allows to control the relaxation process of the excited state of the cluster.

Full size image

The existence of multiple crossing points (A1, A2 and A3) in Figure 2 upon the variation of onsite interaction strength U of both ground and excited states reveals an important implication of the role of frustration in correlated fermionic systems. Not only it leads to a multitude of possible ground-state phases, but also to a rich manifold of metastability patterns upon the decay from excited states. The considered system can be seen as a constituting block for the anisotropic Hubbard model on a triangular lattice, which receives substantial interest being prototypical to the spin-liquid physics of the correlated charge-transfer organic salts60. Richness of the phase diagram is related to an interplay between frustration effects and the anisotropy of the triangular lattice; the latter is present in our consideration as V ≠ V′. We therefore expect rich metastable relaxation dynamics patterns, controlled by the frustration of the lattice in other frustrated strongly correlated systems, which represents a substantial interest for further studies.

The presented approach requires operations with the density matrix and can be directly applied to the correlated cluster of only a small size, due to the exponential growth of the Hilbert space of the system. Certain increase of the system size can be achieved when the truncation of the many-body spectrum is possible. This is the case for systems where a limited number of transition channels works61,62 or at small temperatures of the bath when only the lowest energy states are involved in the dynamics. The relevant part of the many-body spectrum can be efficiently obtained by the Lanczos method. For a construction of the many-body density matrix, there is a variety of established approximate methods of computational quantum chemistry of molecules, including Hartree-Fock method and its extensions63.

To summarize, we considered a correlated triangular cluster where the interplay between interaction and frustration results in level crossings which have an effect on the decay behaviour. We have shown that the dynamics of open correlated quantum systems requires an extensiton of the standard Lindblad one-particle formalism. Working with a many-body spectrum of the correlated triangular cluster, we were able to capture its non-trivial multiscale relaxation dynamics to the thermal equilibrium as a result of weak coupling to the reservoir. We demonstrated the dynamics of the spin and charge degrees of freedom of the system and explained the arising metastability by selection rules which are imposed on transitions between many-body states. The multistage relaxation process involves the excitation to the states lying higher in energy. As a result, the time scale of dynamics becomes dependent on the temperature of the reservoir. The approach can be applied to a number of problems in quantum chemistry, nanoelectronics and quantum information technology.

Methods

In the master equation approach a weak hybridization of an open quantum system with the environment is treated in a perturbative manner. At the same time correlations are taken into account exactly by calculating many-particle eigenstates |n〉 and eigenenergies En of the isolated system. Small reservoir excitations from the equilibrium Gibbs state ρB normally decays fast on a resolved time-scale, so the Markov approximation holds. In the lowest non-vanishing order the resulting equation reads

Substituting HE, HI in this expression is covered in details in Supplementary Materials. For a practical purpose we discard the dependence on energy of the hybridization function assuming that a wide band approximation is fullfilled.

Transition operators are defined as follows

where is a Fermi distribution function. The master equation takes the form

There is a relation of Lα, to the Lindblad formalism approach employing operators and cα in the master equation. The following correspondence was found: 1) , 2) Lα = 0, and 3) Lα = cα, . Though for concreteness the dynamics of fermions is considered in the paper, the developed formalism can be applied to both kinds of particle statistics.