{'fraction_non_alphanumeric': 0.03518552036199095, 'fraction_numerical': 0.006914027149321267, 'mean_word_length': 4.50647797488539, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A complete quasiclassical map for the dynamics of interacting fermions

We present a strategy for mapping the dynamics of a fermionic quantum system to a set of classical dynamical variables. The approach is based on imposing the correspondence relation between the commutator and the Poisson bracket, preserving Heisenberg's equation of motion for one-body operators. In order to accommodate the effect of two-body terms, we further impose quantization on the spin-dependent occupation numbers in the classical equations of motion, with a parameter that is determined self-consistently. Expectation values for observables are taken with respect to an initial quasiclassical distribution that respects the original quantization of the occupation numbers. The proposed classical map becomes complete under the evolution of quadratic Hamiltonians and is extended for all even order observables. We show that the map provides an accurate description of the dynamics for an interacting quantum impurity model in the coulomb blockade regime, at both low and high temperatures. The numerical results are aided by a novel importance sampling scheme that employs a reference system to reduce significantly the sampling effort required to converge the classical calculations.


I. INTRODUCTION
Molecular simulation is an indispensable tool for understanding many-body quantum systems driven away from equilibrium. Describing the dynamics of molecular-and mesoscopic-electronics on time-and length-scales relevant to experiments, however, is challenging. In recent years, significant progress has been made by introducing numerically converged techniques, such as methods that rely on real-time diagrammatic sampling techniques [1][2][3][4][5][6][7][8][9] , wave function-based approaches such as numerical renormalization group techniques [10][11][12] and multi-layer multiconfiguration methods, 13,14 or reduced and hierarchical density matrix approaches. 8,15 While significant progress has been made using these methods to understand the transport in various correlated scenarios, [16][17][18] their application to more realistic systems is still limited.
Previous attempts to map the dynamics of fermionic systems onto a set of classical dynamical variables failed to reliably reproduce correlation effects, such as the Coulomb blockade staircase. 36,[38][39][40] This is mainly due to the lack of quantization of the number operators in the classical map, leading to a continuous increase of the current with the increase of bias or gate voltage, in a quantum point-contact setup. Moreover, the description of the dynamics of observables that depend non-linearly on a pair of creation and annihilation operators, for example, in shot-noise measurements, has not received any attention. As shown below, a naive and straightforward application of the classical maps to such observables leads to significant errors, even for noninteracting model Hamiltonian, where the map generates the exact dynamics.
In this study, we develop a new strategy to map the dynamics of an open quantum system driven away from equilibrium onto a set of classical dynamical variables. The method maps a pair of creation or annihilation fermionic operators to phase-space variables in Cartesian coordinates that satisfies a correspondence relation between the commutator and the Poisson brackets. In order to accommodate the effect of twobody terms (electron-electron interactions), we further impose quantization rules on the spin-dependent occupation numbers in the classical equations of motion, with an onset parameter that is determined self-consistently. Combining this map with the initial value representation [44][45][46] that incorporates the discrete nature of quantum mechanics results in a robust description of the dynamics on diverse time-scales, as illustrated for the Anderson impurity model 47 for a wide range of temperatures and on-site electron-electron repulsion term. We further show that for quadratic Hamiltonians, higher order fermionic operators can be mapped accurately as a consequence of completeness, providing a framework to study the fluctuations and high order correlations within this mapping approach. Finally, we develop a reference sampling approach to reduce significantly the number of trajectories required to converge expectation values.

II. ANDERSON IMPURITY MODEL
For concreteness, throughout this manuscript, we consider the evolution of observables for the Anderson impurity model. This model is defined by the Hamiltonian H = H S + H B + V , whereĤ describes the impurity (or dot), referred to simply as the 'system Hamiltonian',Ĥ describes the noninteracting fermionic baths (or leads), and describes the hybridization between the system and the leads.
Here, d † σ (d σ ) are the creation (annihilation) operators of an electron on the dot with spin σ =↑, ↓ with a one-body energy ε σ . U is the on-site Hubbard interaction,ĉ † kσ (ĉ kσ ) are the creation (annihilation) operators of an electron in mode k of the leads with energy ε k , and t k is the hybridization between the dot and mode k in the lead. The coupling to the quasicontinuous leads is modeled in the wide band limit. The spectral function of the left ( = L) or right ( = R) lead is: where Γ determines the coupling strength to the -lead, B is the width of the spectral function, and A determines the sharpness of the cutoff. The coupling t k between the dot and the k-th mode is expressed in terms of the spectral function as t k∈ = J (ε k ) ∆ε/2π, where ∆ε = 2ε max /(N /2 − 1) is the discretization of the leads energy spectrum, N is the numbers of modes in the -lead, and 2ε max is energy range in the leads.
To model accurately the wide band limit, one should consider sufficiently large values for B such that the energy scale of the system is encompassed inside the spectrum of the leads and that the modes in the leads are dense enough, i.e. ∆ε is sufficiently small. In the simulations below each lead consists of N = 600 modes, where half are with spin up and the other half with spin down. Throughout, we takeh, k B and the charge of the electron e, to be 1.
To assess the accuracy and robustness of the quasiclassical mapping procedure, we focus on the Coulomb blockade effect that is manifested by a staircase structure of the current versus voltage, as shown in Fig. 1. When the bias voltage is not sufficiently large to overcome the on-site repulsion energy, only one conductance channel is open. When the bias becomes sufficiently large compared to U, an additional conducting channel opens up, and the current increases to its maximal value of a two-channel quantum point-contact. In Fig. 1, we show the results of two quasiclassical mapping procedures. The mapping approach that is isomorphic to quaternions (Li-Miller map (LMM)) 38 provides an accurate description of the I-V characteristics at low and high bias voltages (V SD ), but fails to reproduce the Coulomb blockade staircase. On the other hand, the current complete quasiclassical map (CQM) provides a qualitative description across all values of V SD . In particular, it captures the staircase structure characteristic of the Coulomb Blockade effect. We will return to discuss the CQM results after we introduce the strategy of mapping quantum to classical degrees of freedom.

III. COMPLETE QUASICLASSICAL MAP (CQM)
For an operatorÂ, in the Hilbert space of the Anderson impurity model, the Heisenberg equation of motion reads whereĤ 0 , (6) is the one-body, noninteracting part of the Hamiltonian. We wish to find a map forÂ to a function of classical phase-space variables, A[ R], that would preserve the dynamics Â (t) = A(t) c , where Â (t) = Tr ρÂ(t) , and is the classical expectation value with respect to the initial probability distribution ρ 0 of the total system. We do this in two parts. First we construct a complete map for quadratic Hamiltonians, which is extendable to any even order operator, and valid for noninteracting fermionic systems. Then, we propose a strategy for mapping Hamiltonians of higher order containing onsite Hubbard interactions. In all the simulations, the equations of motion are solved numerically using an adaptive Rung-Kutta(4,5) method. The number of trajectories, N tr , used to converge the results is specified below for each case study. For the steady-state results, additional time averaging is considered.

A. Noninteracting fermions
We first consider the case of non-interacting fermions, U = 0, described by a Hamiltonian that depends quadratically on the creation and annihilation operators,Ĥ 0 . To reproduce the dynamics of the expectation value of quadratic operators under the evolution describe in Eq. (5) we require: (a) For any quadratic operatorÂ, and its classical counterpart A, the commutator and the Poisson bracket satisfies the correspondence relation i Ĥ 0 ,Â = {A, H 0 }.
(b) For any quadratic expectation value, the initial probability distribution ρ 0 must satisfy Â (0) = A(0) c , and respect the quantum discrete nature of the occupations.
It is straightforward to show that condition (a) is satisfied by mapping a pair of creation and annihilation operators to a phase space of conjugated variables, R = (x, p x , y, p y ), aŝ a † nân → x n p y,n − y n p x,n , ,n x m + y n p y,m − p y,n y m ) + (x n p y,m − p x,n y m + x m p y,n − p x,m y n )] , ,n x m − y n p y,m + p y,n y m ) − (x n p y,m − p x,n y m − x m p y,n + p x,m y n )] , This identifies positions, (x, y) and their conjugate momenta, (p x , p y ), and the Poisson bracket can be used to check that this map returns the quantum commutator of any pair of quadratic creation/annihilation operators. Because any quadratic Hamiltonian with a set of quadratic operators constitutes a closed Lie algebra of quadratic operators, condition (a) insures a loyal representation of the dynamics in terms of Hamilton's equations. We note in passing that we subtracted 1/2 from the classical map ofn i ≡â † iâ i to include a Langerlike correction 48 .
For leads that are in thermal equilibrium and uncorrelated initial state, condition (b) can be satisfied by setting the initial occupation of each mode in the left and right leads to a value 0 or 1, such that the expectation value, averaged over the set of initial conditions, satisfies the Fermi-Dirac distribution. 36 Operationally, we choose a random number ξ kσ ∈ [0, 1] and then select the occupation of mode kσ of the -lead according to where β = 1/k B T and µ l are the inverse temperature times Boltzmann's constant and chemical potential of the -lead, respectively. The Cartesian coordinates are then sampled according to 39 where θ kσ is chosen randomly in the interval [0, 2π] and n kσ = x kσ p y,kσ − y kσ p x,kσ satisfies Eq. (9), resulting in n kσ = n kσ c = (1 + exp[β (ε k − µ )]) −1 at the initial time. By construction, the expectation value at the initial time is a † n a m c = a † n a † m c = a n a m c = 0, as expected for averages taken with respect to uncorrelated thermal distribution. The sampling choice in Eq. (10) is not unique, but it does provide an efficient averaging of the expectation values with respect to the number of trajectories. 39 In a similar manner, one can set the initial occupation of the dot.
Comparing the proposed CQM given by Eq. (8) to the LMM, we find that the mappings of the diagonal termâ † nân and of the linear combinationâ † nâ † m +â † mâ † n are identical in both maps, but the remaining terms in Eq. (8) cannot be expressed using the LMM. This leads to the Hamiltonian being expressed identically in both maps, For this non-interacting Hamiltonian, mapping H 0 and then deriving Hamilton's equations of motion for the phase space variables is identical to deriving Heisnberg's equation of motion for the bi-linear operators and then mapping the results using Eq. (8).
We can also map other quadratic observables, such as the current from the left lead: ∑ σ =↑,↓ k∈L t k y σ p y,kσ − p y,σ y kσ + x σ p x,kσ − p x,σ x kσ .
As a diagonal term, the above form is also identical to the expression obtained by the LMM. 39 In the upper panel of Fig. 2 we compare the results for the left current generated by the CQM (which in this case are equivalent to the LMM) with exact quantum mechanical results for a noninteracting model Hamiltonian. As expected, the agreement between the CQM (or the LMM) and exact quantum mechanical results is excellent. In the next section we show that for the CQM these result can be extended to higher-order operators.

B. Higher order operators
Mapping higher order operators, operators that involve more than one pair of creation/annihilation operators, is more difficult due to the nonlocal character of fermions arising from their exclusion statistics. Ignoring the fermionic nature does not seem to make any significant difference for a single pair of creation/annihilation operators, 36,39 but for higher order operators, the anti-commutation of the creation/annihilation fermionic operators plays a significant role and describing quantum fluctuations such as shot noise requires a careful consideration of this effect.
For example, consider a map for the operatorÂ = a † n a m a † m a k . Using the anti-commutation nature of {a † n , a n } = 1, we can express the expectation value ofÂ using four different terms that are identical quantum mechanically, but differ when mapped onto classical phase space variables. Specifically, expanding Â we find there are four unique combinations of operators, which generically have coefficients, C i . To determine the best choice of C i , we impose conditions (a) and (b) of Sec. III A on the time evolution ofÂ and require that the time evolution ofÂ be exact for a quadratic Hamiltonian, i.e., that i[Ĥ,Â] = {A, H} and that Â (0) = A(0) c . For an uncorrelated initial thermal state, the values that satisfy these conditions are C 1 = 1,C 2 = −1,C 3 = 1, and C 4 = 0. Note that Eq. (13) contains pairs of creation or annihilation operators (a † n a † m or a n a m ), which cannot be described within the LMM.
Applying this procedure to the second moment of the left current for the noninteracting Hamiltonian yields a simple ex- As can be shown explicitly, this mapping of the second moment of the left current operator satisfies both conditions (a) and (b).
In the bottom panel of Fig. 2  between the exact quantum mechanical result and the CQM is excellent. We also plot the individual terms Î 2 L,1 , Î 2 L,2 , and Î 2 L,3 (dashed lines). Only the proper combination of all three terms yields an accurate description of Î 2 L . We note that the LMM can only be used to map the first term, but not the other two that contribute to Î 2 L .

C. Interacting fermions
The on-site interaction, Un ↑n↓ , that manifest the Coulomb blockade effect, is a four index term that is outside the space defined by the CQM. In order to account for this two-body interaction term, we map the two terms proportional to U in Eq. (5) according to where θ is the Heaviside step function. The idea behind this choice is that the term Un ↑n↓ contributes to the dynamics only when both electrons with spin up and spin down occupy the site. Classically, the occupation number admits a continuous value, which implies that the Hubbard term can become significant for fractional populations of the two spin-channels. Much like a mean-field approximation, these fraction contributions of the Hubbard term will smear the Coulomb blockade effect. By introducing the step function, we impose that contributions to the dynamics from the Hubbard term arise only in trajectories for which n ↑(↓) > ∆ ↑(↓) . The parameter ∆ ↑(↓) is determined according to the distribution of n ↑(↓)c and will be discussed in detail below. We note that the classical expression in Eq. (15) is not derivable from a Hamiltonian, and therefore does not in principle conserve energy or the norm of phase space. Nevertheless, we find relaxation to an intermediate time, long lived steady-state for all of the observables studied on timescales shorter than the recurrence times.
Considering the Anderson impurity model, the equations of motion for the Cartesian variables for the lead degrees of freedom areẋ jσ = −ε j y jσ − t j y σ , (16) y jσ = ε j x jσ + t j x σ , p x, jσ = −ε j p y, jσ − t j p y,σ , p y, jσ = ε j p x, jσ + t j p x,σ , and those for the system's degrees of freedom arė whereσ =↓, ↑ is the opposite spin to σ =↑, ↓.
In Fig. 1 we plot the I-V curve obtained by the CQM and compare it to the QME in Ref. 49 and to the results obtained by the LMM. We consider the limit of weak system-bath coupling and high temperature where the QME provides a good approximation for the dynamics of the system. The LMM provides a good description of the I-V characteristics at low and high bias voltages, but it fails to capture the staircase structure reminiscent of the Coulomb blockade. The CQM reproduce the QME results quantitatively, specifically, it captures the staircase structure due to the Coulomb blockade effect. In this high temperature regime, the agreement between the CQM and the quantum mechanical results is observed for a wide range of onsite Hubbard repulsion term and also for the quantum dot population.
Next, we consider a regime where the QME breaks down, namely, the low temperature regime. For simplicity we focus on the equilibrium case where µ L = µ R = 0. In this regime, solutions for the population as function of the gate voltage (ε = ε ↑ = ε ↓ ) are readily available using the numerical renormalization group (NRG) technique. 50,51 In the left panel of Fig. 3 we plot the quantum dot total population ( n ↑ + n ↓ ) as a function of the gate voltage. The NRG results show a staircase shape which is a manifestation of the Coulomb blockade effect. The CQM agrees quantitatively with the NRG results over a wide range of gate voltages. Specifically, it captures both the position of the blockade as well as its width. The QME approach, however, captures only the position of the resonances; the broadening of the transitions are missing completely. This qualitative difference between the CQM and QME approaches signifies the advantages of the quasiclassical mapping techniques over the commonly used QME approach for a broad range of temperatures.
The mapping of the Hubbard term in Eq. (15) introduces a parameter, ∆ ↑(↓) , which is determined self-consistently. For the results show in Fig. 3 we use a single value for ∆ σ = 0.18, determined by considering the particle-hole symmetry point, where ε ↑ = ε ↓ = −U/2. At the symmetric point, the steady state value of the average dot populations is n σ = 1/2. Since quantum mechanicallyn σ can only assume two value, 0 or 1, whereas the distribution of n σ is continuous, to obtain an average dot population n σ = 1/2, we set ∆ σ to the median of the distribution of values of n σ . This ensures that the Hubbard terms in Eq. (15) are significant only when n σ is sufficiently large, modeling the discrete nature of the spin-dependent dot occupations.
In the right panel of Fig. 3 we show the results obtained for the total dot occupation for different values of ∆ σ . Only 3 iterations are required to converge the results for ∆ σ . We start with an initial guess of ∆ σ = 0.34. In the next iteration, we set ∆ σ to the new median of the distribution of values of n σ , in this case ∆ σ = 0.24. We then repeat this procedure until convergence (only one more iteration is required). The results clearly show that the expectation value is not very sensitive to small variations in the value of ∆ σ , but the converged results provide the best agreement with the NRG results.
The agreement between the CQM and the quantum mechanical results are not limited to steady-state properties. In fact, our mapping also captures quantitatively the hallmarks of the Coulomb blockade in the relaxation towards the steady state. Shown in Fig. 4 is the time dependence of the dot occupation for two different values of U, with a bias voltage of V SD = µ L − µ R = Γ, and temperature T = Γ. Here we compare the CQM to numerically converged real-time stochastically sampled diagrammatic techniques applied within the reduced density matrix formalism. 7 For both values of U we find that the full time dependence is in good agreement with the numerically converged data. In each case, the dot population increases monotonically and the timecales required to reach the long time limit are comparable.

IV. REFERENCE DYNAMICS FOR STATISTICAL CONVERGENCE
When evaluating the classical expected value Eq. (7), the integral over the initial distribution is replaced by averaging over different initial configurations of the leads that satisfy the Fermi-Dirac distribution. For a large number of initial conditions, N tr , the procedure converges to the desired distribution and to the exact expectation value. However, the low dimensional nearly harmonic system generically requires a large number of initial conditions to statistically converge the result. To reduce the number of initial conditions for a given statistical error, we introduce a reference system whose expectation value can be determined exactly and inexpensively. Specifically, the expectation value of an observable A is calculated according to where A r is an observable used as the reference, andĀ r is the exact expectation value of A r evaluated using a different inexpensive method. In the limit N tr → ∞ we have A r →Ā r and A will approach the real expectation value. However, for a finite N tr , and a smart choice of A r , one can reduce significantly the statistical error of this estimator. This is clarified by considering the variance using the reference system If A and A r are correlated, it is possible to have Var(∆A) < Var(A). As we are now propagating both A and A r , to reduce the computational effort we desire that Given that the sample variance reduces as 1/N tr , to obtain computational superiority the ratio in Eq.(20) is bounded by 1/2. However, improvement in the computational effort can already be seen for ratios that are above this factor, since typically propagation of the reference system is not as costly as of the system of interest. Non-interacting or mean-field Hamiltonians that can be solved analytically serve as examples of reference systems that can reduce noise for the dynamics of interacting systems. Other possibilities include considering dynamics that are generated from some effective Hamiltonian with the same initial configurations. Fig. 5 are two examples, one in which the reference system method works good and one in which it fails. On the left panel we plot the ratio of the variances of the left current with spin up as a function of time, and on the right we plot the corresponding currents. The reference considered here is the current calculated for noninteracting systems where an exact solution can be obtained by direct diagonalization of the single-particle Hamiltonian. We see that when the reference current,Ī r , becomes very different from the real current, I r + ∆I , (for U = 6Γ in the figure) the ratio of the variances at steady-state exceeds its bound 1/2. However, when the reference and the real currents are proximate but still quite different (for U = Γ/2 in the figure), we see that the ratio of the variances is reduced significantly as a consequence of correlations between the trajectories of the currents.

Shown in
We note that for parameter regime where the reference and real currents almost coincide the fluctuations drastically decreased, the ratio of the variances at steady state reaches as low as ∼ 10 −4 . This implies that for a fixed statistical convergence threshold, the number of initial conditions decreases by two orders of magnitude, since each estimate is statistically independent. One can also note that at short times the reference system always reduces the fluctuations significantly. The reason is that we used an uncorrelated initial condition and thus the short time behavior is set by ∼ Γ −1 . It takes a certain amount of time for correlations to build up and for the interacting part in the Hamiltonian to influence the dynam-ics, yielding a statistical benefit for short times even when the steady-state result is far from the noninteracting limit.
The idea of using a reference system can be extended beyond the description above. For example, if one wishes to calculate the current as function of U, one can start the evaluation for small U and increase it "adiabatically". For each calculation of the current, the previous current (with smaller U) can be used as the reference system. Of course, the exact term in Eq. (18) is no longer exact and carry with it some error, but this can still be beneficial, as trajectories of the different currents are likely to be correlated given that the change in U is small.

V. CONCLUSIONS
We have presented a quasiclassical method to simulate nonequilibrium dynamics of interacting fermions. We have constructed this map using the correspondence relation between the commutator and the Poisson bracket, in order to preserve Heisenberg's equation of motion for one-body operators. We have shown that this classical map is complete for quadratic expectation values under quadratic Hamiltonians and it can be extended to higher moments accurately. This feature makes the study of fluctuations and higher-order correlations accessible.
For interacting systems, the dynamics is approximated by mapping the equation of motion and enforcing a quantization rule that determines for which values of n σ the dynamics is influenced by the Hubbard term. This, together with a quasiclassical initial distribution, provides a quantitative agreement with other methods in regimes where those other methods are known to be accurate. Thus, a quantitative description of nonequilibrium currents in the Anderson model, including their steady-state behavior as illustrated by the presence of the Coulomb blockade, their fluctuations as encodes in the second moment of the current, and the relaxation of each to their steady state can be obtained.
We have also shown a way to enhance the statistical convergence of this method by introducing a reference system, whose dynamics can be computed exactly, and averaging the difference between the reference system and the system of interest. Provided the reference system is correlated with the system of interest, fluctuations are reduced in the averaging procedure, and we have shown that this can increase the computational efficiency by up to 2 orders of magnitude over naive sampling. Together, these results make the quasiclassical method appealing for studying nonequilibrium phenomena in complex chemical systems. Indeed, realistic systems of molecular junctions routinely operate a low effective temperatures, and finite interaction strengths rendering other low scaling approximate approaches inaccurate. The method we have presented here is capable of probing these regimes, at a small computational cost that scales linearly in the system degrees of freedom. This should enable studies in correlated transport behavior in high dimensional, molecular systems, far from equilibrium. Using the relations in Eqs. (21) and (22), quadratic creation and annihilation operators can be expressed as: The commutation relation of two elementary quaternions are then mapped to a cross product of vectors in phase space, The CQM replaces the map in Eq. (24) with  24) can be used to map the termsâ † nân andâ † nâm +â † mân . The CQM extend this to terms likeâ † nâ † m ,â nâm andâ † nâm .