Many-body Tunneling and Nonequilibrium Dynamics of Doublons in Strongly Correlated Quantum Dots

Quantum tunneling dominates coherent transport at low temperatures in many systems of great interest. In this work we report a many–body tunneling (MBT), by nonperturbatively solving the Anderson multi-impurity model, and identify it a fundamental tunneling process on top of the well–acknowledged sequential tunneling and cotunneling. We show that the MBT involves the dynamics of doublons in strongly correlated systems. Proportional to the numbers of dynamical doublons, the MBT can dominate the off–resonant transport in the strongly correlated regime. A T 3/2–dependence of the MBT current on temperature is uncovered and can be identified as a fingerprint of the MBT in experiments. We also prove that the MBT can support the coherent long–range tunneling of doublons, which is well consistent with recent experiments on ultracold atoms. As a fundamental physical process, the MBT is expected to play important roles in general quantum systems.

Although the doublon have been studied in ultracold atom systems, the connection to the MBT and its many-body nature are not fully appreciated. We demonstrate that the doublon dynamics can be described by two-particle density matrix and is a kind of MBT naturally originating from strong Coulomb interaction. We select multi-impurity Anderson models 23 for illustration due to the tunable system parameters, e.g. double, triple and quadruple quantum dots (DQDs, TQDs and QQDs) connected to two biased reservoirs. The models are nonperturbatively solved by the hierarchical equations of motion approach (HEOM). The intrinsic features of the MBT and doublon are studied and their universality is addressed.

MBT in DQDs
Fundamentally, both ST and CT are only related to the reduced time-dependent single-particle density matrix, ρ(t), thus their current is not sensitive to U. By contrast, the MBT engages reduced two-particle density matrix, π(t). There are altogether three kinds of π(t), corresponding to three kinds of basic physical processes (see the equations of the HEOM in Ref. Jin08234703). They are: 1) the Kondo  (doublon) and its Hermitian conjugation (holon). In these expressions, c is ( + c is ), with s = ↑, ↓, denotes the annihilation (creation) operator matrix of an electron in the specified spin state in the i th dot; = + n c c is s s is the number operator with opposite spin of s; a js and b js denote dot or bath operator matrixes which can be the same or different. In the present work, neither the Kondo effect nor the doublon-holon excitation plays a role within the investigated parameter range, thus we precisely work out the MBT current.
Although it is accompanied by the production and decay of doublons, based on the analytical analysis, it is noted that the MBT only concerns the tunneling of electrons, rather than the direct tunneling of doublon-holon excitations. The latter is another important many-body effect which will be addressed in future works. We also comment that all of the physical quantities obtained from HEOM calculations, including current, occupation number and spectral functions, will be self-consistently checked to achieve their precise values in what follows. Our numerical results prove that no other contributions to the current besides ST, CT and MBT, which confirms that the MBT involves all high-order tunneling processes beyond the second order.
By referring Fig. 1(a), one can deduce that the Kondo spin screening may occur if spin flips during the doublon decay 24,25 . In single QD, the Kondo resonance dominates the transport in the impurity magnetic moment regime, overwhelming the MBT by competition. However, the MBT may make leading contribution in multi-QD systems. First of all, we introduce Pauli spin blockade (PSB) in DQDs. The PSB results from the Pauli's exclusion principle which prevents the spin triplet from transporting out of the DQDs, and is characterized by the electron density accumulation as well as the current suppression in one bias direction [26][27][28] . The reason of choosing the PSB regime is that the MBT dominates the transport leading to a finite current even when the CT is blocked. In fact, the current can be large enough to degrade the PSB effect, which can be a measurable symbol of the MBT ( Fig. 2(a)).
To make the discussion concrete, we numerically calculate the nonequilibrium spectral functions A is (ω) in the PSB regime, as shown in Fig. 1(b) and (c). The MBT involves four events occurring almost simultaneously. i) A ↓-spin electron tunnels from the left reservoir to QD1; ii) The excess ↓-spin tunnels from the QD1 to QD2; iii) The ↓-spin electron in the QD2 excites to the doubly occupied level (doublon production); and iv) The doubly occupied ↓-spin electron in the QD2 tunnels to the right reservoir (doublon decay). A key point is that the MBT strongly depends on U in both doublon production and decay processes. In comparison, the CT shown in Fig. 1c) can be well described by the second-order time-nonlocal quantum master equation (QME) 29 . Different to the MBT, the tunneling from the QD1 to the right reservoir is now via the singly occupied level of the QD2, without doublon production nor decay process. This leads to insensitivity of the CT to U, which is an essential difference of the CT to the MBT.
The Anderson multi-impurity model is adopted for N-QD (N = 2, 3 or 4) systems. The total Hamiltonian reads H total = H sys + H res + H sys−res . The Hamiltonians of reservoirs are c ks ) denotes the annihilation (creation) operator of an electron in the specified spin state in the α-reservoir with wave vector k. We set The hybridization function is assumed to be a Lorentzian form , with W = 4 meV fixed in this letter. The QD-reservoir coupling strength Γ will be specified below. The Hamiltonian for central N-QD is For DQD case, we lift the spin degeneracy of the QD1. This can be achieved with a local spin-splitting micromagnet 30  . For the tilted TQDs (see the insert of Fig. 5(a)), E 1 = −E 3 = E, E 2 = 0; and for the tilted QQDs (see the insert of Fig. 5(c)), The multi-impurity Anderson model is solved using the HEOM approach 31,32 . This approach supports accurate and efficient evaluations of various steady-state and transient properties [32][33][34][35] (for more details, please refer to Ref. [YeWIREs] and references therein). It has been demonstrated that the HEOM approach achieves the same level of accuracy as the latest high-level numerical renormalization group and quantum Monte Carlo approaches for the prediction of various dynamical properties at equilibrium and nonequilibrium 31 .
The current-voltage characteristics at t = 0.01 and 0.01 meV are shown in Fig. 2(a) and (b), respectively. For comparison, the QME results essentially from the perturbative treatment are shown. Let us examine closely the forward current, I(V > 0). In both t = 0.01 meV and t = 0.01 meV cases, the QME currents show negative differential conductances, due to the transition from the resonant to off-resonant tunneling 36,37 . At large V, the QME currents tend to near-zero constant value resulting from the CT, with a ratio of I peak /I const ≈ 50. For precise HEOM currents (including the nonperturbative MBT currents), things are different. In t = 0.01 meV case, the MBT current is increased and the negative differential conductance almost vanishes. Although in the t = 0.1 meV case, both the QME and HEOM currents display similar negative differential conductance, it should be pointed out that the HEOM current in Fig. 2(b) gives us I peak /I const ≈ 4, which matches the experimental data well 36 . Although some other external mechanisms concerning the leakage current in the PSB regime have also been proposed in literatures 38 , here we provide a more intrinsic and universal one. In contrast to I CT ∝ t 2 resulting from the second-order perturbation theory, I M exhibits a strong nonlinearity before saturation. By fitting the numerical data of I M in Fig. 2(d), we find I M ∝ exp[−α(U)/t], where α(U) is a parameter approximately proportional to U. It means that the MBT current has a similar dependence on U/t to the decay rate of doublons observed in Ref.
[Strohmaier10080401], which originates from the many-body character of the doublon dynamics. It is more important that these features are rather general, beyond the case of eV = U in Fig. 2

(c) and (d).
For spin-nondegenerate-QD1 case, the tunneling current is basically carried by down spins, i.e., ↓  I I , with ↑  I 0. As shown in Fig. 1, ↓ I M should be closely related to the number of doublons generated in QD2 and decayed into the right reservoir. We call this part of doublons as 'dynamical doublons' and denote their number as n DD . In the DQDs, n DD (t) = 1/2 − n 2↑ (t), where n 2↑ is the ↑-spin occupation number in QD2. For the uncoupled limit n 2↑ (t = 0) = 1/2 leading to n DD (t = 0) = 0. ↓ I M is evaluated as a function of n DD for different Γ and U in Fig. 3 (a) and (b), respectively. Remarkably, ↓ I M keeps a fundamental linear function of n DD , i.e.  Fig. 2(c). I M ∝ n DD proves that the dynamical doublons are the only carriers of I M , namely, the MBT precisely describes the dynamics of doublons, as we argued above.
The ideal linear plots in Fig. 3(c) shows further scaling laws of the slope k on the ratio of Γ/U, k ∝ (Γ/U) 2 . In Fig. 3(d), we prove the temperature dependence of the MBT current as, 3/2 , which is distinctly different from the T 2 -dependence of I CT 9 . This unconventional T dependence, originally attributed to the intrinsic many-body character of the MBT, can be conveniently verified by experiments.
It is necessary to extend the study to the spin-degenerate case in which the local Zeeman splitting in QD1 is removed. We thus choose = = ↓ ↑ 0 1 1   and keep other parameters unchanged. The corresponding I-V curves are shown in Fig. 4(a), the counterparts to Fig. 2(a). It is observed that the MBT current still dominates in the weak coupling regime at high bias. Apparently, the main features of the MBT current in the spin-degenerate case are basically the same as those of the spin-nondegenerate case. Figure 4(b) is the spin-degenerate counterpart to Fig. 2(c). Again the basic features remain unaltered. In particular, I M also increases exponentially with t 2 , and then gradually saturates at a constant value. From the above results, we can conclude that the MBT characters still hold well in the spin-degenerate case. In other words, the MBT is a general process.

Long-range MBT
Now, we are on the position to demonstrate the coherent long-range MBT in larger systems than DQDs. In recent experiments of the Bose-Hubbard chain at t ≪ U, long-range tunneling over up to five sites were observed as resonances in the number of doublons when the nearest-neighbour tilt energy (E) is tuned to integer fractions of U 7 . Let us prove that the observations arise from the long-range MBT by numerical results in tilted TQDs and QQDs (see inserts in Fig. 5(a) and (c)). Figure 5(a) depicts the dependence of I, I M and I CT on t 4 at eV = 1.5 meV for the spin degenerate TQD with E = U/2. Not surprisingly, all of them behave in a similar manner as those in DQD. Since the number of tunneling barriers doubles now, we have I CT ∝ t 4 .
In tilted QDs, I CT can not be blocked as that in the PSB region, thus it is mixed with I M in observed I. In order to separate out the direct contribution of the MBT, we define γ ≡ (I M /I × 100%). Figure 5(b) shows the dependence of γ on the ratio of E/U for the TQDs. The case of titled DQDs is shown in the insert for the purpose of comparison. For the TQDs, a resonant peak at E = U/2 is clearly seen; while for the DQDs, that peak shifts to E = U, with the peak position unchanged at different t. The maximum value of γ up to 75% in TQDs further proves that the resonance at E = U/2 results from the long-range MBT rather than the CT which in fact has been excluded already by the strong U-dependence of the resonance. Both the peak structure and the peak position for the DQDs and TQDs shown in Fig. 5 Our theoretical results for titled QQDs further confirm the above arguments. In Fig. 5(c), I CT is approximately proportional to t 6 at E = U/3. In Fig. 5(d), a similar peak structure is clearly seen with the resonant peak located at E = U/3, which is exactly the position of resonant doublons tunneling through four sites of titled Hubbard chain, as observed in expriments 7 . We thus conclude that the MBT can support the coherent long-range tunneling of doublons in general systems, which is of both fundamental and practical importance. For example, the MBT can provide a feasible way to manipulate distant quantum gates or qubits in one step in solid-state quantum computing. This process would enhance the operating efficiency and fault-tolerant capability, compared to the nearest-neighbor control in exchange-based quantum gates 39 .

Summary
In summary, we have theoretically discussed the many-body tunneling (MBT), by nonperturbatively solving the Anderson multi-impurity model, and identified it a fundamental tunneling process that involves the dynamics of doublons. Proportional to the numbers of dynamical doublons, the MBT is shown to dominate the off-resonant transport in strongly correlated systems. A T 3/2 -dependence of the MBT current is uncovered and can be identified as a fingerprint of MBT in experiments. It is also proved that the MBT can support the coherent long-range tunneling of doublons. Our theoretical results are well consistent with recent experiments on the dynamics of    , where Δ is the effective quantum dot-electrode coupling strength, W is the band width, and μ α is the chemical potentials of the α (α = L, R) lead.
In this paper, a hierarchical equations of motion approach (HEOM) developed in recent years is employed to study QD system 31,32 . The HEOM based numerical approach is potentially useful for addressing the interacting strong correlation systems and has been employed to study dynamic properties, such as the dynamic Coulomb blockade Kondo, dynamic Kondo memory phenomena and time-dependent transport with Kondo resonance in QDs systems. The resulting hierarchical equations of motion formalism are in principle exact and applicable to arbitrary electronic systems, including Coulomb interactions, under the influence of arbitrary time-dependent applied bias voltage and external fields. The outstanding issue of characterizing both equilibrium and nonequilibrium properties of a general open quantum system are referred to in ref. 40. It is essential to adopt appropriate truncated level to close the coupled equations. The numerical results are considered to be quantitatively accurate with increasing truncated level and converge.