QuTiP-BoFiN: A bosonic and fermionic numerical hierarchical-equations-of-motion library with applications in light-harvesting, quantum control, and single-molecule electronics

The"hierarchical equations of motion"(HEOM) method is a powerful exact numerical approach to solve the dynamics and find the steady-state of a quantum system coupled to a non-Markovian and non-perturbative environment. Originally developed in the context of physical chemistry, it has also been extended and applied to problems in solid-state physics, optics, single-molecule electronics, and biological physics. Here we present a numerical library in Python, integrated with the powerful QuTiP platform, which implements the HEOM for both bosonic and fermionic environments. We demonstrate its utility with a series of examples. For the bosonic case, we include demonstrations of fitting arbitrary spectral densities, and an example of the dynamics of energy transfer in the Fenna-Matthews-Olson photosynthetic complex, showing how a suitable non-Markovian environment can protect against pure dephasing. We also demonstrate how the HEOM can be used to benchmark different strategies for dynamical decoupling of a spin from its environment, and show that the Uhrig pulse-spacing scheme is less optimal than equally spaced pulses when the environment's spectral density is very broad. For the fermionic case, we present an integrable single-impurity example, used as a benchmark of the code, and a more complex example of an impurity strongly coupled to a single vibronic mode, with applications to single-molecule electronics.


I. INTRODUCTION
In the context of open quantum systems [1] many approaches exist to modelling the influence of an environment on a system.For example, when the interaction between system and environment is weak and Markovian (memory-less), a range of perturbative methods are available, primarily in the form of generalized Lindblad master equations [2,3], or non-secular master equations [1,4].However, beyond these perturbative limits one must start to treat the dynamics of the environment on the same level as the system, which in many cases is a challenging task.
Discretization of a continuum environment is one powerful approach to this non-perturbative limit [5][6][7][8][9][10], wherein the predominant collective degrees of freedom are identified and included in the full simulation in a numerically efficient way.The "hierarchical equations of motion" (HEOM) method does precisely this, albeit indirectly [11][12][13][14][15].It is based on the construction of a hierarchy of coupled equations resulting from taking repeated time-derivatives of an influence functional, under the assumption that the bath correlation functions take an exponential form.Such an exponential form can be acquired analytically for certain spectral densities (with Matsubara or Padé decompositions [16,17]) or numerically with fitting [18].Once the set of equations is obtained, they can be numerically solved to give the system dynamics, the system steady-state, and certain bath properties.The limitations of the HEOM arise in the truncation of the set of equations (which in the case of a bosonic environment is infinite) and the truncation of the exponential decomposition of the bath correlation functions.The latter limits the accuracy to represent the real problem and can lead to unphysical results if poor decompositions are made.
In this paper we introduce an open-source library for modelling these equations of motion for both bosonic and fermionic baths as an integrated part of our larger toolbox for modelling open quantum systems in Python, QuTiP [19][20][21][22].We describe the bosonic (Sec.II) and fermionic (Sec.III) forms separately.For both we start with basic definitions, and then provide a range of examples, some of which reproduce important known results, and three of which give novel insights gained with this tool for this article: • In Sec.II C we apply the HEOM to the ubiquitous spin-boson model.We show results for the Drude-Lorentz spectral density, and the under-damped Brownian motion spectral density, and show convergence trends for the Matsubara, Padé, and fitting approaches to decomposing the bath correlation function.We also compare the results to a standard Markovian method in the form of the Bloch-Redfield solver in QuTiP.
• In Sec.II D we demonstrate how to model multiple baths by reproducing the seminal results of [13] for the dynamics of energy transport through the Fenna-Matthews-Olson pigment protein complex, and showing that the non-Markovian nature of the environment preserves electronic coherence as compared to a standard Bloch-Redfield approach.
We also take the recent experimental results of [23] and, using a combination of the HEOM and the reaction-coordinate method, show that systemenvironment entanglement oscillations persist on longer time-scales than electronic coherence in a strongly-damped regime.
• In Sec.II E we illustrate how the auxiliary density operators in the HEOM can be used to obtain information about the environment by reproducing results from [24] regarding non-equilibrium heat flow into an environment.
• In Sec.II F we demonstrate how capturing the non-Markovian nature of the environment can be important for quantum control of quantum systems by examining the interplay between dynamical decoupling schemes and environment properties.Importantly, we both demonstrate how the HEOM can be used in regimes where standard analytical results fail (finite control pulse length), and benchmark two pulse-spacing schemes against each other (equally spaced versus Uhrig's optimal strategy).
• In Sec.II G we show how to go beyond the standard choices of spectral density with the HEOM using different fitting techniques.We explicitly demonstrate how to capture an Ohmic spectral density with exponential cut-off using either spectraldensity fitting or correlation function fitting techniques, and benchmark them against each other using the pure-dephasing model.
• In Sec.III D we benchmark the fermionic solver against an analytical result for the current through a single resonant level coupled to two fermionic baths.
• In Sec.III E, for the fermionic solver, we reproduce the results in [25] showing the steady state current through a resonant level coupled both to two fermionic baths and to a single vibronic mode.
Each example is associated with a Jupyter Notebook provided on GitHub [26], which throughout this paper are referred to via the numbered index in their name, 1a, 1b, etc.

A. Comparison to existing packages
Other open-source implementations of the HEOM exist, a summary of which, including the one presented in this work, is given in table I A. It is difficult to be exhaustive on the capabilities and what constitutes useful levels of documentation or testing for all these packages.Hence, for documentation, a "yes" was given if the available documentation seemed sufficient for a third-party to install and run the software on their own problems, without having to read the source code.A "yes" was given for tests if the software included a suite of tests that cover the main features of the implementation and checks the resulting outputs.QuTiP's implementation is the only one that meets this standard.
In summary, Tanimura provided an early Fortran implementation on his website [27].PHI [28,29] was a later implementation that provided multithreaded CPU support.The Nanohub GPU HEOM [30] (unrelated to the Tsuchimoto GPU-HEOM [31]) provided GPU support but is closed source.The Tsuchimoto GPU-HEOM [31] (also from the Tanimura group) implements a custom matrix exponential solver to reduce memory overhead.HEOM-QUICK [32] is the only other implementation to support fermionic baths, and provides much flexibility via editing functions in the Fortran code.DM-HEOM [33][34][35] supports distribution across multiple nodes, allowing the simulation of systems where the memory requirements of the hierarchy exceed that available on a single node.PyHEOM [36] provides explicit support for a new decomposition method for capturing poles in the bath spectral densities.
Our implementation is complementary to the above packages, in the sense that it provides flexible implementations of both bosonic and fermionic cases, and allows for both pre-defined commonly used spectral densities and arbitrary user-based input (that can be combined with custom fitting of bath correlation functions or spectra, as required).In addition, integration with the QuTiP library will ensure its continuous maintenance and improvement by an established team of developers [37].

A. Basic Definitions
When considering the influence of an environment on a quantum system, a typical starting point is to consider the environment as a bath of linear harmonic oscillators.In some physical situations this is justified by the actual dominant degrees of freedom being harmonic in nature.However, it is also a very pragmatic and often phenomenological assumption, as, when combined with the assumption that the bath is initially in a Gaussian (e.g., thermal) state, evaluating the influence of the bath on the system is much easier.In this work we do not consider non-linear baths (e.g., spin-baths), though it has been shown that such environments can be captured with the HEOM to some degree [38,39].
In the standard second-quantized Hamiltonian formalism, we can write the interaction of a single bosonic environment with an arbitrary system in the following way, where H S (t) is the free Hamiltonian of the system, which can be time-dependent, H B = k ω k a † k a k is the free Hamiltonian of the bath, Q is the system operator which couples to the bath, and g k is the coupling strength between the system and mode k in the bath.Note that this, and the following, is easily generalized to multiple baths with different system coupling operators.For convenience, we define Note that in the above, and throughout most of this work, we set = k B = 1 (some later examples are presented in different units).The state of the system at time t has an exact solution in the form of a time-ordered influence functional, which only depends on the free bath correlation functions [1,11,12,40], Here ρS (t) is the density matrix of the system in the interaction picture with respect to H S (t)+H B (indicated by the bar) at time t, and there is an implicit assumption that at t = 0 the system and the environment are in a product state ρ(t = 0) = ρ S (t = 0) ⊗ ρ B (the sub-index B refers to the bath or environment).The environment is an initially thermal equilibrium state with temperature T = 1/β and Z = Tr[e −βH B ]. Importantly, the system operators Q that couple to the environment act to the right as superoperators in Eq. ( 3), with The derivation of this exact result, Eq. ( 3), relies crucially on the Gaussian nature of the free environment, which allows the properties of the system to only depend on the second order correlation functions of the environment.These correlation functions can be expressed as, where n(ω, β) = 1/[exp(βω) − 1], and in moving from the discrete to continuum limit one defines Returning to Eq. ( 3), which is formally exact but still a time-ordered exponential that is difficult to directly solve, if we assume that an appropriate choice of J(ω) gives correlation functions with real and imaginary parts that can be decomposed as, where c j k and γ j k themselves can be real or complex, one can formally take repeated time derivatives of Eq. ( 3) to arrive at an infinite set of coupled first-order differential equations.A concise description of this procedure can be found in the appendix of [40] for a particular case (see also [41] for an alternative approach, and [11] for the original derivation using path integrals), but for the most general decomposition one ultimately finds the following set of coupled differential equations [42], where where we use the same notation for the commutator as in Eq. (5).In some cases this can be augmented with additional Lindblad terms, as described later.
If there are terms in the decomposition of the correlation function with equal frequencies γ R k = γ I k they can be combined into a single index, for a gain in numerical efficiency [42].This is done automatically by our code itself.Due to the linearity of the interaction with the bath, additional environments can be included by simply extending the list of bath correlation functions (i.e., adding a bath label K to Eq. ( 9), and different system coupling operators to different baths is allowed for by labelling the indices and the coupling operator Q K appropriately.

B. Code Functionality
Our Python-based library relies on, and is integrated with, the popular Quantum Toolbox in Python (QuTiP) [19,20], and hence uses the ubiquitous Qobj data types from that library to represent states, operators and superoperators.Generally speaking, one sets up the problem to be solved by defining the system Hamiltonian (which can be time-dependent), the bath properties using either pre-defined bath class objects, or through lists defining the correlation function decomposition given in Eq. ( 9) and a generic bath class, and a system coupling operator for each bath.
The HEOMSolver class constructs the right-hand-side of Eq. (11), taking as input one or more baths (which can be defined in multiple ways), the system Hamiltonian, and the truncation parameters.For example, for a bosonic bath using the predefined Drude-Lorentz spectral density with N k Matsubara terms: bath = D r u d e L o r e n t z B a t h (Q , lam = lam , gamma = gamma , T =T , Nk = Nk , tag = " DLBath " ) HEOM_dlbath = HEOMSolver ( Hsys , bath , NC )

Code Listing 1. Defining the hierarchy of equations
Here Q is the system-bath coupling operator, and the other parameters determine the properties of the Drude-Lorentz spectral density Eq. (15).In addition, one can define a bosonic bath using the same spectral density with a Pade decomposition [DrudeLorentzPadeBath()], the underdamped spectral density of Eq. ( 16) [UnderDampedBath()], or a generic one using the real and imaginary coefficients of Eq. (39) [BosonicBath()].
This HEOMSolver object can then be used to solve the coupled set of ordinary differential equations using standard libraries, or calculate the steady-state.For the former purpose, the class contains a run() function which takes an initial condition and a set of time-steps, and returns a results object with both the system density matrix and the auxiliary-density-operators (ADOs) for the bath at each time-step: result_dlbath = HEOM_dlbath .run ( rho0 , tlist , ado_return = True )

Code Listing 2. Solving the time evolution
The state of the system is returned for each time step in result_dlbath.states,while, if ado_return=True is used, then the full state of the auxiliary density operators is returned as a list as well, in result_dlbath.ado_states.These can then be used to obtain certain bath properties.Tags can be employed to delineate different baths, and levels can be used to extract different tiers.For example, for the example defined above, ado_states.filter(level=1,tags=["bathDL"]) returns a list of labels identifying the ADOs that are tagged with "bathDL" and are at truncation tier 1 (i.e., only have j=R,I Nj k=1 n jk = 1 for the given bath).Using result_dlbath.extract(),these labels can then be used to extract specific ADOs, see Example 3 and its associated notebook for a practical example.
For computing the steady-state solution, we provide two direct methods, one of which takes advantage of Intel's MKL library if installed alongside QuTiP.Full details of the functionality can be seen in the accompanying Jupyter notebooks (provided in [26]) and documentation [43].

C. Example 1: Spin-boson model
The archetypical test system for studying open quantum systems is a two-level system (i.e., a spin, or qubit, but we use the terminology two-level system throughtout this work).In the context of physical chemistry, it is used as a model of electronic energy transport between two nearby molecules (i.e., a dimer model).Referring back to Eq. (1), in this example we describe the system Hamiltonian using the standard Pauli operators to represent the two-level system and the coupling operator Ultimately the choice of spectral density, J(ω), depends on the properties of the environment one is considering.A useful choice is that of an Ohmic (linear) frequency dependence with a Lorentzian cut-off.This is usually split into two types, the over-damped Drude-Lorentz form and the under-damped Brownian motion form, While superficially similar, the properties of the corresponding correlation functions (C D (t) and C U (t)) are different in several important ways.Firstly, as we will show below, with a Matsubara decomposition C D (t) has only overdamped exponents, while C U (t) can have oscillatory parts.Secondly, J D has a 1/ω high-frequency tail, while J U converges much more quickly.This leads to two unique features in C D (t) at t = 0: a non-zero imaginary part (which is arguably unphysical [44]), and a divergent real part.For example, the Matsubara decomposition of the Drude-Lorentz spectral density is given by: The divergent real part at t = 0 can be captured in the HEOM formalism by treating exponentials above a certain cut-off, N k , as delta functions [12].Since It is convenient to calculate the whole sum and subtract off the contribution from the finite N k Matsubara terms that are kept in the hierarchy, and treat the residual as a prefactor multiplying a delta function.It is then possible to show [12] that the delta-function contribution to the correlation functions can be described by Lindbladian terms in the HEOM; and thus the equation of motion in Eq. ( 11) is modified so that the Liouvillian evolution includes these Lindblad terms where This treatment is sometimes called the "Tanimuraterminator", and an example of how to include it is presented in the accompanying example notebook 1a[26].
As described earlier, the QuTiP-BoFiN package [45, 46] will automatically determine that real and imaginary exponents are close, and if the coupling operator for both exponents are equal, it will combine them into a single common exponent.Note that for degenerate exponents in the real or imaginary terms alone the code will not automatically combine terms.
1. Two-level system coupled to a Drude-Lorentz bath In Fig. 1 we show a simple example of the dynamics of a two-level system coupled to a Drude-Lorentz bath.We compare several optional approaches.This includes an explicit truncation of the Matsubara decomposition at N k = 2, the same truncation with the Tanimura terminator, and an example where a large number of Matsubara terms (N k = 15×10 3 ) are summed up and numerically fit using the fitting function C F = N f i=1 a i e bit with N f = 4 exponents.We also show the solution from the standard Bloch-Redfield solver in QuTiP (i.e., a generalized Born-Markov-Secular master equation).In the accompanying example notebook 1a [26] we also present the results of a Padé decomposition of the correlation functions, but do not show it in the figure (it essentially converges faster than the Matsubara case alone as it gives the same result as the Matsubara plus terminator decomposition, but without the need of an equivalent terminator).
We can see from Fig. 1 that, for this choice of parameters, either the Tanimura-terminator, the Padé decomposition, or the fitting approach is needed to reach the correct steady-state (alternatively more discrete terms must be included in the Matsubara summation).The black dashed lines show the steady-state obtained using a reaction-coordinate approach [7].This trend can be confirmed with a pure-dephasing analytical result, which is provided in example notebook 1e [26], but not explicitly shown here.
2. Two-level system very strongly coupled to a Drude-Lorentz bath Moving to a more challenging parameter regime (see also example notebook 1b [26]), Fig. 2 shows the situation of a very strong coupling between system and environment (as used as a benchmark in [17]).We compare the results obtained with a single Matsubara term, the same with the terminator included, a single Padé decomposition term, and a case where we sum up a large number of Matsubara terms and fit them with auxiliary exponents.Here the Bloch-Redfield approach fails completely, so it is not explicitly shown.As described in [17], in this case the terminator result converges very quickly, while many discrete Matsubara terms would be needed without it.The Padé result also converges faster than the single-term Matsubara result, as does the result using a fit, but with a higher numerical overhead.We expect a more sophisticated fitting algorithm may provide better results (as evidenced by the better performance of the Padé decomposition with less exponents).

Two-level system coupled to an underdamped bath
Moving to the underdamped case, the underdamped spectral density gives correlation functions (again using the Matsubara decomposition) separated explicitly into real and imaginary parts of the form: Where Ω = Ω 2 c − (Γ/2) 2 .Note that here we slightly abuse notation to include the pair of complex conjugate terms with one index (k = 0).These are presented as separate terms in the real and imaginary decomposition here, but the package will again combine real and imaginary parts of the k = 0 components with equal exponents, such that there are two total effective exponents instead of four.
In Fig. 3 (see also example notebook 1c [26]) we show a typical result for a strongly coupled, narrow spectral Figure 1.HEOM results for the excited state probability ρ11 and the coherence ρ01 versus time for a two-level system coupled to a bosonic environment described by a Drude-Lorentz spectral density.The system parameter used here is = 0.5∆, and the bath parameters are λ = 0.1∆, γ = 0.5∆, and T = 0.5∆.We compare a solution of the Bloch-Redfield master equation to the Matsubara decomposition with N k = 2, the same decomposition with the Tanimura terminator, and the case where only the primary (Drude) exponent (k = 0) is kept in Eq. ( 17), and a large number (N k = 15 × 10 3 ) of Matsubara terms are summed and fit with N f = 4 exponents.The horizontal dashed black line shows the thermal (steady) state result obtained from a reaction-coordinate approach.The code for generating this figure can be found in example notebook 1a [26].
density.We see a large deviation between the Born-Markov-Secular Bloch-Redfield result and the HEOM result, and a convergence of the HEOM result with around four Matsubara terms.In particular, we observe that the Matsubara terms, for narrow spectral densities like this one, largely contribute to correcting detailed balance and arriving at the correct steady state.We benchmark this effect by plotting the steady-state using the reaction coordinate (RC) method [7], which shows how as we add more Matsubara terms to the HEOM approach we tend towards the steady-state predicted by the RC method (see [8] for a more detailed analysis of this effect, and results that include fitting of Matsubara terms to reach T = 0).Matsubara HEOM results for the excited state probability versus time for a two-level system coupled to a bosonic environment with a Drude-Lorentz spectral density in the strong coupling regime.Here we use the system parameter = 0.0, ∆ = 0.2γ, and the bath parameters λ = 2.5γ, and T = γ.We show the case of a single Matsubara exponent, the same with the Tanimura terminator, the Padé decomposition also with a single (non-Drude) exponent, and again a fit of N k = 15 × 10 3 Matsubara terms, this time with N f = 3 exponents.The code for generating this figure can be found in example notebook 1b [26].

D. Example 2: FMO-complex
The HEOM method has a long history of being applied to problems in physical chemistry, particularly in the study of transport of electronic excitations through light-harvesting complexes [13,[47][48][49].The typical parameters of such systems sit in a regime where it is hard to justify any particular perturbative approximation, as the electronic dipole coupling between chromophores and the dominant environmental influence are, in many cases, of the same order.One of the most studied examples of this type of problem is that of the Fenna-Mathews-Olson (FMO) complex.It serves us well here as a way to show how to use our package for a system with several internal levels coupling to independent environments.
One again starts with a Hamiltonian description of the system.We model a single excitation moving through the FMO complex with a seven-site model where the on-site energy of site i is given by i and the dipole coupling between bacteriochlorophyls (BChl) i and j is given by J i,j .Again we treat the environment as a continuum of harmonic oscillators, and assume that each BChl is coupled to an independent but identical bath of such oscillators.This represents how the harmonic nuclear motion of the molecules leads to changes in the electronic energies.The coupling operator for the i-th bath is thus Thermal state Figure 3. HEOM result for the excited state probability versus time for a two-level system coupled to a under-damped spectral density with system parameter = 0.5∆, and bath parameters α = 0.5∆ 3/2 , Γ = 0.5∆, ωc = ∆ and T = 0.05∆.Here we show how the result converges as we increase the number of Matsubara terms, and compare to the solution from the Bloch-Redfield master equation solver (which fails drastically).The black dashed line represents the thermal steady-state predicted by the reaction coordinate method.The code for generating this figure can be found in example notebook 1c [26].
described by a Drude-Lorentz spectral density.A separate bath class object can be created for each bath, associated with each system coupling operator, using whatever bath-decomposition scheme one requires.An explicit example of this is given in example notebook 2 [26].
We provide an example of simulation results in Fig. 4 showing how, with the initial condition of an excitation localised on site 1, the populations evolve up to 1 ps.These results identically recreate those seen in [13].
To show how a non-Markovian environment differs from a Markovian one, in Fig. 5 we show the results, using the same parameters, for the standard Bloch-Redfield solver in QuTiP.This solver, under standard Born-Markov-Secular approximations, takes into account the frequency dependence of the bath spectral density, the eigenstate structure of the system Hamiltonian, and the pure dephasing that arises from the low-frequency behavior of the bath.
In the accompanying notebook 2 [26] we explicitly show how these different contributions affect the dynamics, and we see there quite explicitly that the pure dephasing (in the eigenbasis) has a very strong effect on the dynamics, largely suppressing coherent oscillations.This comparison explicitly shows that weak-coupling models fail to predict the correct coherence timescale in this case, due to the highly non-Markovian environment [13,50]  Figure 4. HEOM results for the site-populations of the FMO complex versus time.The FMO complex consists of seven sites coupled to seven identical but independent baths each described by a Drude Lorentz spectral density with parameters λ = 35 cm −1 , γ −1 = 166 fs, T = 300 K.For the FMO Hamiltonian we use the data employed in [13].In comparison to Fig. 5 we see that the exact solution predicts much longer coherent oscillations than the equivalent weak-coupling solution.The code for generating this figure can be found in example notebook 2 [26].
determined partially by the narrow cut-off γ).In earlier works suppression of low-frequency noise was shown to arise for super-Ohmic spectral densities [51], based on arguments regarding the frequency dependence of the Markovian pure-dephasing contribution, but appears here due to the non-Markovian nature of the environment.Interestingly, whether the electronic coherence, as exhibited in the oscillations seen in Fig. 4, actually survives in nature is still hotly debated [23].In [23] the authors observe no electronic coherence beyond 60 fs, in contrast with earlier experiments [52][53][54][55].In modelling their data they employ a theoretical method using a much broader and stronger bath spectral density than that used in [13] that strongly suppresses the coherent oscillations of the type shown in in Fig. 4.Here we adapt the spectral density employed in the supplementary information of [23]: with the parameters listed there [56] (see Fig. 6(a) and (b) for a plot of this function).
To gain some insight into this parameter regime, we first truncate the FMO space to only the first two-sites coupled to a single bath.This gives a reduced total Hamiltonian similar to our earlier examples, with a coupling operator to a single bath Q = σ z / √ 2, where the factor √ 2 arises in the reduction of the original two-bath model to the normal modes of a single-bath model (see [57] for more details on this step in the reduction).Using the FMO parameters we employed earlier we find = −120 cm −1 and ∆ = −87.7 cm −1 , and again use T = 300 K.
We then simulate this reduced system using the HEOM in two ways: first we do a very rough approximation using a single overdamped Drude-Lorentz spectral density, then we do a more nuanced fit using two underdamped spectral densities.We then repeat the simulations using the reaction-coordinate (RC) approach [7], a method that is more approximate than the HEOM but known to perform well at high temperatures [58] even with broad spectral densities.
Importantly, the RC method gives us access to a density operator for the original system plus an effective collective "reaction coordinate" mode, from which we can obtain system-environment properties less directly accessible with the HEOM.In addition, since the RC approach amounts to performing a Bogoliubov transformation on the original bath modes combined with tracing out a residual environment, we can calculate the system-RC entanglement via quantities like the negativity.We expect this to be a lower-bound on the entanglement between the system and the original full environment (since said entanglement cannot be increased by local operations and classical communication), if the perturbative approximation for deriving a master equation for the residual bath modes holds.
Figure 6 demonstrates the results from these various approaches.First of all, we use the exact HEOM method to validate the more approximate RC method in this parameter regime, and we see that for both choices of spectral densities the methods produce identical system population dynamics and steady states, and that both spectral density choices produce essentially the same result.This appears to be because the large coupling at low frequencies dominates the system dynamics and causes quick suppression of coherent dynamics in both the site [ρ 11 in Fig. 6(c) and (d)] and eigenbases [ρ ge in Fig. 6(c)  and (d)].
Furthermore, Fig. 6(e) shows the result of calculating the negativity between the system and the single RC mode used for the overdamped Drude-Lorentz spectral density, and the two modes used for the twounderdamped spectral densities.We observe dynamics in this quantity up to several 100 fs and significant nonzero steadystate values in both cases.Note that the non-secular master equation used to describe the residual bath of the reaction-coordinates can sometimes induce non-positive eigenvalues in the RC modes themselves, as it does not necessarily preserve positivity.This potentially signals the start of a breakdown of the perturbative approximation needed to derive this master equation, and thus suggests a limit on how much we can infer about the original system-bath entanglement.However, for the parameters used in Fig. 6, negative eigenvalues only appear at very short times and are small compared to the negativity values shown in Fig. 6(e).
In a more complex model involving more sites, we expect the negativity between the first two sites and the environment to decay with time as the electronic excitation moves through the complex.Nevertheless these results suggest that observing electronic coherence alone can be insufficient to determine whether non-trivial systemenvironment correlations can be neglected, particularly when the effective coupling to the environment is on the order of the system energies.In addition, we see here the benefit of combining the exact HEOM with more approximate but nuanced methods, like the RC approach.

E. Example 3: Quantum heat transport
In the examples discussed so far, we have investigated the time evolution of the system density matrix alone.An intriguing feature of the HEOM is that it also allows us to study observables that depend on the state of the environment, which is encoded in the auxiliary density operators (ADOs).In this example, we will demonstrate this feature by calculating the quantum heat current in a two two-level system setup.The heat transport in this setup was originally studied by Kato and Tanimura in Ref. [24]; we here show that their results can be reproduced with ease using our framework.
The setup consists of two coupled two-level systems, which are each in turn coupled to an individual heat reservoir.The Hamiltonian describing the two-level sys-   (b) shows a fit instead using two under-damped brownian motion spectral densities, giving α 2 1 = 0.11ω 3 c,1 , ωc,1 = 900 cm −1 , Γ1 = 800 cm −1 , and α 2 2 = 0.25ω 3 c,2 , ωc,2 = 500 cm −1 , Γ2 = 990 cm −1 .(c) shows the HEOM and RC results obtained for dynamics of several matrix elements on the two-site dimer model using the fit from (a), while (d) shows the same using the fit from (b).The RC and HEOM produce comparable results, and the choice of fitting function has little influence on the system dynamics in this case.(e) shows the entanglement between system and RC mode(s) as determined by the negativity using the single RC mode model for the fit in (a) and and two RC modes using the fit in (b).
tems is where denotes the level splitting, J 12 the coupling between two-level systems, and σ K z,± denotes the usual Pauli matrices for the K-th two-level system.Each reservoir is modeled as a continuum of harmonic modes with an overdamped Drude-Lorentz spectral density as defined in Eq. ( 15).The multi-reservoir setup can be easily treated J 12 = 0.01 J 12 = 0.1 J 12 = 0.5 Figure 7. HEOM results for the steady-state heat current in the two two-level system setup.The parameters were chosen as follows: γ1 = γ2 = 2 , T1 = 2.02 , T2 = 1.98 and λ1 = λ2 = ζ(J12/2), where the dimensionless parameter ζ parametrizes the system-bath coupling strength.This plot is a reproduction of the HEOM-curves in Figs.3(a-i) to 3(a-iii) of Ref. [24].The code for generating this figure can be found in example notebook 3 [26] .
using the HEOM by including a bath index K ∈ {1, 2} in the multi-index of the ADOs.The bath coupling operators are given by Q K = σ K x .A temperature difference T 1 > T 2 induces heat transport in the setup.Due to the second law of thermodynamics, the heat flow is on average directed from the hot side to the cold one.Our goal is to determine the individual heat currents between the system and the reservoirs, which are generally given by j being the Hamiltonian of the K-th reservoir.These can be determined from the ADOs as follows [59]: We here added indices signifying the reservoirdependency to the rates γ j;K k , the cut-offs N K j , the correlation function C K (t) and the Tanimura-terminator Γ K T .The symbol ρ (jk;K) denotes the level-1 ADO corresponding to the multi-index (0 . . .0, 1, 0 . . .0) with only a single non-zero entry at the indicated position.
In our library, the HierarchyADOsState class facilitates the extraction of specific ADOs from the full state of the simulation.The formula in Eq. ( 32) can therefore be easily implemented, as we demonstrate in detail in example notebook 3 [26].In Fig. 7, we show the behavior of the steady-state heat current j ss = j 2,ss = −j 1,ss in the hightemperature regime as a function of the system-reservoir coupling strength.The figure shows that the heat current vanishes both for large and small coupling strengths and reaches a maximum in the intermediate regime.The behavior at large coupling strengths is caused by a quantum Zeno-like effect, where a strong system-reservoir interaction suppresses the coherence of the system, inhibiting the energy flow between the individual two-level systems [24,60,61].

F. Example 4: Driven systems and dynamical decoupling
Another powerful aspect of the HEOM lies in the fact that it is valid for any system Hamiltonian, even timedependent ones.As a concrete example, we show here how it can be used to study dynamical-decoupling, a common tool used to "undo" dephasing from an environment, even for finite pulse duration (i.e., away from the "bang-bang" pulse limit).
The general idea of dynamical decoupling is that dephasing from a non-Markovian environment can be undone by the application of certain choices of fast pulsecontrol on the system itself.We follow the protocol outlined in [62], which involves rapid applications of π pulses, interspersed with periods of interaction with the environment.The intuitive concept is that if these pulses act faster than the memory time of the bath, by swapping the state of the two-level system, the effect of the environment can be reversed due to a change in sign in the interaction.
Here we choose again a Hamiltonian in the form of Eq. ( 13), but setting ∆ = 0 and adding a time-dependent drive term, We use the same coupling operator as Eq. ( 14), and the same spectral density as Eq.(15).We move to an interaction picture with respect to the time-independent part of H S , such that the drive Hamiltonian has the form, The pulse is chosen to be, for equal spacing ∆ t between pulses, V n (t) = V , for n∆ t + (n − 1)τ p ≤ t ≤ n∆ t + nτ p , and zero elsewhere, (see Fig. 8b for a schematic of the pulse shapes).The period of the actual pulses themselves is τ p V = π/2.Note that in comparison to [62] we omit some phase factors that appear if one considers a more realistic model of the drive.In Fig. 8a (see also example notebook 4 [26]) we show the time-evolution of the coherence of the two-level system assuming the initial condition ψ(t = 0) = 1 √ 2 (|0 + |1 ).When the τ p is chosen to be fast (and the corresponding V strong so that V τ p = π/2) we see the expected cancellation of dephasing occurring (green curve).If the pulse is too slow, we see the onset of dephasing (blue curve).The decoherence expected with no pulse at all is also shown (orange dashed curve).
It has been shown [63,64] that for environments with very sharp cut-offs [65] (e.g., step-functions), a squaredsinusoidal choice of spacing is optimal (a scheme referred to as Uhrig Dynamical Decoupling (UDD)).In that case, one sets V n (t) = V when and zero elsewhere.This produces pulses that are clustered more closely together at the beginning and end of the overall time-evolution period T max .The Drude-Lorentz spectral density we use for this example has a very long tail, and thus this approach is not guaranteed to be optimal at all times.In Fig. 9 we demonstrate this by comparing the final coherence ρ 01 at time t = T max as a function of γ and for 100 decoupling pulses.We see that for smaller γ the UDD outperforms the equally spaced case, but as we increase γ its performance drops.Interestingly, as expected from the discussions in [63,64], the coupling strength λ does not reduce the relative effectiveness of the UDD appreciably.These results demonstrate how HEOM can be used to validate control schemes in realistic scenarios such as finite pulse length and realistic choices of spectral density.
To evaluate time-dependence in the Hamiltonian, our code mirrors the typical approach taken for standard solvers used in the QuTiP package.The time-dependence is defined through a function that is sent alongside the dependent part of the Hamiltonian in a list.For the above example, this follows simply as

H = [Hsys, [sigmax(), drive]]
where drive is a function determining the timedependence.The various parameters can also be defined as arguments for the function.Our code also accepts QobjEvo format for the time dependence.

G. Arbitrary spectral densities
A great deal of effort has been made in the literature to optimize the analytical decomposition of the Drude-Lorentz spectral density into exponentials, so as to minimize the numerical overhead in the HEOM.For more general spectral densities, one can consider two approaches: fitting the spectral density directly with Drude-Lorentz or underdamped spectral densities, or evaluating the correlation functions and fitting those directly.Both approaches have been discussed in the literature [66,67], with the general consensus being that the optimal choice depends on what one is evaluating -the dynamics of short or long time scales, the spectral response of the system, etc.A more recent method also suggests generalizing the types of functions used in the decomposition of the correlation functions [36].
To compare the different fitting approaches, we can consider a class of spectral densities with polynomial fre- For the slow drive we choose V = 0.02 Vf , and in both cases we choose a pause period of ∆t Vf = 10.The lower figure, (b), shows the drive amplitudes as a function of time.In (a) we see that very slow pulses perform much worse than fast pulses, as expected due to the competition with the bath memory time.Note that here and in Fig. 9 we operate in a rotating frame with the drive, and assume that the system frequency is much larger than the drive amplitude ( V ) such that the rotating wave approximation is valid.The code for generating this figure can be found in example notebook 4 [26].
quency dependence and exponential cut-off: The correlation functions for these spectral densities are [68] wherein Γ is the Gamma function and ζ is the generalized Zeta function.
Here we present one simple example of fitting the Ohmic case s = 1, but the following can also be applied to other examples.If one chooses to fit Eq. ( 36 .Example of dynamical decoupling showing the coherence ρ01, after a total evolution time Tmax = 1000(1/ V ), of a two-level system coupled to an environment with a Drude-Lorentz spectral density versus γ, the width (inverse memory time) of the bath.Each run includes 100 dynamical decoupling pulses, and we compare equally spaced pulses (labelled DD in the figure legend) with the UDD scheme described in the text.The parameters here are chosen as ∆ = 0, λ = 0.001 V , 0.0005 V , 0.00005 V , and T = 0.05 V .
There is a clear crossover for larger values of γ where the UDD scheme becomes less efficient than the standard equally spaced scheme (DD).The code for generating this figure can be found in example notebook 4 [26].
directly, Meier and Tannor [18] proposed using a set of underdamped spectral densities in the form Here, in converting this decomposition into the same form as those in Section II C, we set Γ i = Γ i /2.In Fig. 10 (see also example notebook 1d [26]) we show an example of such a fitting, for an Ohmic environment s = 1, using k J = 4 terms.With the fitting parameters we can also reconstruct the correlation functions and power spectrum, as a check of validity of the fit.
Instead, if we choose to fit the correlation functions, we treat the real and imaginary terms separately, and expand them as In Fig. 11 we show the fitting for the same parameters as Fig. 10 using k R = k I = 3 terms, as well as the reconstructed power spectrum.One can see that at negative frequencies this fit has some residual oscillations, which may lead to some small error in the dynamics.In both spectrum and correlation function fitting, it is important to restrain the fitting parameters to some degree, as arbi-trary parameters can make achieving convergence in the resulting HEOM simulation difficult.
We compare the results of both approaches in Fig. 12 for parameters close to Fig. 2, except with slightly narrower bath and longer time-scales.More importantly, for the spectral density decomposition, we employ a terminator to compensate for the relatively low truncation of the Matsubara terms (see Fig. 10).The number of ADOs, and hence numerical cost, with the correlation fitting method is 2(k R + k I ) .For the spectral density fitting each k J term adds two exponents for the non-Matsubara term, and a single exponent for each matsubara term.Hence the number of exponents is k J (2 + N K ).In Fig. 12 we see that k R = k I = 3 gives comparable results to k J = 3, N k = 1 (N k = 0 performs badly and is not shown), implying similar numerical cost for approximately the same accuracy in this case (with slightly less exponents needed for the spectral fitting approach).However, as one lowers the temperature, more and more Matsubara terms are needed in the spectrum approach.This is corroborated by Fig. 13 which compares the fitting parameters in the HEOM against the analytical result available when we set the system Hamiltonian to commute with the bath coupling operator (here achieved by setting ∆ = 0).First we see that, in Fig. 13(a), for the same parameters and temperature as Fig. 12, for k J = 3 and k R = k I = 1 we see similar magnitude of accuracy, while k R = k I = 3 does exceptionally well.In Fig. 13(b), at a temperature lowered by a factor of a half compared to Fig. 12, we find that the spectral decomposition requires both more spectral terms and more Matsubara terms to be on par with the correlation function fitting result.In particular, the cyan curve for k R = k I = 3 requires twelve exponents in total, while the blue curve for K J = 4 and N K = 4 requires 24 exponents, suggesting much lower numerical cost for correlation function fitting in the low temperature regime.
In all cases we use a standard least-squares algorithm for the fitting procedure.More sophisticated fitting approaches [67], or taking advantage of known properties of the bath correlation functions, may produce better overall results.One may also consider fitting the power spectrum as a third method (which equates to fitting the Fourier transform of the correlation functions).Also, while one can use analytical limits in some case to benchmark the results (like the pure dephasing result in Fig. 13), in general one may lack a way validate any given result.One potential strategy is to perform both fitting procedures (spectrum and correlation functions), and check convergence of the system dynamics in both cases (large differences between them suggest one or both contain non-negligible errors).A more formal refinement of this approach may be to develop error bars on a given result based on the different fitting procedures used, akin to the error bars used in [8] based on noise introduced into the fitting procedure.for a two-level system coupled to an Ohmic environment with exponential cut-off [Eq.(36) with s = 1] for α = 3.25, T = 0.5ωc, = 0, ∆ = 0.2ωc.We show results for fitting the spectral density with kJ = 1, and 3 underdamped Brownian motion terms, and for fitting the correlation functions with kr = kI = 1 and 3 terms.For the case where we fit the spectral density, we use N k = 1 and the standard terminator to approximate the other Matsubara terms.Here, and in Fig.

III. FERMIONIC ENVIRONMENTS A. Basic Definitions
In the previous section we summarized how a bath of bosonic modes can influence a discrete system.Another common scenario is that of a macroscopic conductor coupled to a microscopic impurity, in which case one must consider how discrete states interact with a continuum of fermionic modes.Such a is traditionally analyzed with a range of many-body techniques, but in recent years it was shown that the HEOM method can be applied to this scenario as well.Here we summarize the basic definitions, and present some examples from the fields of quantum transport and single-molecule electronics.
Again, in the standard second-quantized Hamiltonian formalism, we can write the interaction of a fermionic mode c, {c, c † } = 1, with support on the system space, with a fermonic bath, where, as with the bosonic case, H S is the free Hamiltonian of the system (which can contain arbitrary degrees .Difference in the dynamics as evaluated by the HEOM and that given by an analytical result for the coherence |ρ10(t) − ρ ana 10 | for a two-level system initially prepared with maximal coherence, and coupled to an Ohmic environment with exponential cut-off [Eq.(36) with s = 1] for α = 3.25, = 0, and ∆ = 0.The latter means that, unlike Fig. 12, this result is integrable, and we can compare to an analytical result.(a) is for T = 0.5ωc, as in Fig. 12, while (b) is for a lower temperature of T = 0.25ωc.In (a) we show results for fitting the spectral density with kJ = 1, and 3 underdamped Brownian motion terms, and for fitting the correlation functions with kr = kI = 1 and 3 terms.For the case where we fit the spectral density, we here use N k = 1 and the standard terminator to approximate the other Matsubara terms.In (b) the lower temperature means we have to both increase the number of spectral functions and employ more Matsubara terms to reach results comparable to the correlation function fitting results.
of freedom) and is the free Hamiltonian of the electronic reservoir.
As with the bosonic case, we can characterize the free environment through its second-order correlation functions.However, now we discriminate between the order in which excitations are created or destroyed in the environment, so that we define two different correlation functions, defined by the choice of σ = ±: In the following we primarily follow the notation and choices made by Schinabeck et al. [25].Unlike in the bosonic case, we do not divide the correlation functions explicitly into real and imaginary parts, but in the form where η and γ can be complex numbers.The hierarchical equations of motion for such an environment [25,69,70] are then written as where now j = (σ, l) is a multi-index, L = −iH × S is the normal system evolution, and as before ρ (0) is the density operator of the system and n > 0 are auxiliary densityoperators (ADOs) .The superoperators are Here we use σ = −σ, c + = c † , and c − = c.The above definitions are easily extended to multiple baths by extension of the index j, as in the bosonic case.Still following the notation of [25] in this generalization, we will refer to the multi-index j = (K, σ, l), such that there are N j = 2N K (l max +1) indices where N K are the number of reservoirs, and the factor of two arises for the two signs of σ, and (l max + 1) are the number of exponents in the decomposition of the reservoir correlation functions.The primary differences between the bosonic and fermionic HEOM are the following.Firstly, in the fermionic HEOM multi-indices j can only take the values 0 or 1, due to the fermionic nature of the operators in the environment, and thus the truncation parameter N C is the maximum amount of non-zero indices in j, and N C ≤ N j ( equality implies no truncation).Secondly, because fermionic operators anti-commute, the ordering of the indices is now important.Hence the various parity terms that arise in Eq. (43) depending on how many nonzero indices are in the list.This is, in particular, vital for the second term in Eq. ( 43), where on raising an index value from 0 to 1, one can notice that we append the new non-zero index to the left of the list, ρ (n+1) jjn...j1 .Thus one must then move the index to its appropriate position in the unique multi-index j.In doing so, we pick up an additional parity on this term depending on how many non-zero terms appear between the position on the left and default position of said ADO in the list.Note that in other implementations [32], an additional optimization is done by taking advantage of the hermiticity relation between ADOs with conjugate σ.We do not implement that optimization at this stage.
To demonstrate the application of this method, we employ a Lorentzian spectral density of the form where µ is the chemical potential, W is the width of the environment, as before, and η a coupling strength.For fermions, the distribution in Eq. ( 41) is We again have the choice of decomposing the correlation functions into exponentials in several different ways.
The Matsubara decomposition proceeds similarly to the bosonic case.Thus, here we present the Padé decomposition, which approximates the Fermi distribution as where k l and l have to be evaluated numerically, and depend on the choice of l max (see [16], and the example notebooks 5a and 5b [26], for details).Performing this decomposition, and evaluating the integral for the correlation functions gives (see supp.info. of [25]) The Padé decomposition we used for the bosonic case follows in a similar way.

B. Observables
As with the bosonic method, one can directly solve for the dynamics and steady-state properties of system observables.Generally in transport problems, for which the fermionic method is particularly useful, other quantities of interest are the steady-state current, conductance, higher-order transport statistics, and spectral properties.The current, for example, depends on properties of the environment, and can be extracted from the auxiliary density operators, rather than just the system density matrix.It is related to the first-order ADOs (see section 2.2.4 of [71] for a detailed derivation): As with the heat-engine examples discussed earlier, the HierarchyADOsState class allows the extraction of specific ADOs from the results.These can then be straightforwardly used with the above formula (see example notebooks 5a and 5b [26]).

C. Code Functionality
The functionality of the fermionic solver is largely the same as the bosonic one.Each bath is specified via its decomposition into correlation functions, and the associated system operator coupling operator should be provided as a single operator associated with the coupling to the d k modes in the bath (i.e., the single fermion c in the examples below).We provide a Matsubara (LorentzianBath) and Pade decomposition (LorentzianPadeBath) of the Lorentzian spectral density given in Eq. ( 46), or a generic bath (FermionicBath).Examples of this functionality can be found in the documentation and the example notebooks [26].

D. Example 1: Integrable single-impurity model
To benchmark the accuracy of the code it is useful to consider an integrable example of the single-impurity Anderson (SIAM) model.In this case we consider a single spin-less fermionic impurity coupled to two reservoirs We assume the reservoirs are described by Eq. (46) and Eq. ( 47) but generalize them to have different chemical potentials so that the impurity sees a bias ∆µ = µ L −µ R .The steady-state current is a well-known analytical result [72], where the Lamb shifts are In Fig. 14 we show the analytical current and the result from the HEOM, using the Padé decomposition, as a function of the bias voltage ∆µ, for η = 0.01 eV, T = 300 K, and W = 1 eV, with the impurity energy = 0.3 eV.From this figure we see quite clearly that the current flow through the impurity is zero until the bias window is equal or larger than .The current then follows quite closely the energy dependence of the reservoir power spectrum multiplied by the Fermi function.We note that, as discussed in-depth in [73], the HEOM here converges at n max = 2.

E. Example 2: Vibronically assisted transport
Moving away from this simple, integrable, case, we consider the example of a single impurity coupled to a vibronic mode.Such a model is widely discussed in the literature of single-molecule electronics and nanosystems [74,75].Here we explicitly reproduce an example from figure 1(a) in [25].The model used there considers Eq. ( 54) with the addition of an explicit vibronic degree of freedom in terms of a single bosonic mode: In solving this example one has several options: treat the bosonic mode as a HEOM bosonic environment (as done in [74]), include it explicitly in the system Hamiltonian with a discrete Fock space representation, perform a polaron transformation first (as done in [25]), or finally, perform a polaron transform and absorb the vibronic coupling into the lead couplings (which is then treated with an additional approximation, as done in [76] and [75]).We choose the second option here, and truncate the bosonic Fock space at a level N which gives convergence.It is equivalent to the first and third options, while being the most general approach for the demonstration of our library, and is also numerically more difficult than the polaron and direct HEOM approach.In Fig. 15 we show an example of the current versus the bias voltage for = 0.3, eV, Ω = 0.2 eV, λ = 0.12 eV, T = 300 K, W = 10 4 eV, and η = 0.01 eV.We see that at large bias voltages we need a large number of bosonic Fock states to see convergence in the current, consistent with the appendix of [25].The result appears to converge to the result reported in that work for l max = 5 and N = 34 and N C = 2.Note that, if one instead employs the polaron transform approach, a higher-tier  eV, λ = 0.12 eV, T = 300 K, W = 10 4 eV, and η = 0.01 eV.Note that increasing lmax tends to increase the current at large bias voltages, while increasing the other convergence parameter, the bosonic cut-off N , tends to decrease it, gradually giving a converged result in the vicinity of the red curve.
The code for generating this figure can be found in example notebook 5b [26].
result (e.g., N C = 3 here) can be obtained at little numerical cost because its affect on the tier below can be found analytically.We do not employ this optimization here, however, as it is not generally applicable to all problem Hamiltonians.Interestingly, while the current seems to converge, the correlation functions obtained with the Padé decomposition are not converged (due to the wideband limit of W = 10 4 eV).This "quicker" convergence of system properties occurs because, for fermionic environments, certain highly over-damped contributions to the bath correlation functions do not contribute to the system dynamics (see the supplementary information of [77,78] for a detailed proof).

IV. CONCLUSIONS
The HEOM method has, in the last ten years, grown in popularity and applicability.The package we present in this work [43,45,46] captures most of this general applicability in being able to treat both bosonic and fermionic environments, the modelling of multiple environments and time-dependant system Hamiltonians, as we have demonstrated above.We have reproduced seminal results presented elsewhere, particularly the FMO treatment of [13] and the vibronic transport of [25], and also demonstrated some novel applications, including dynamical decoupling with finite length pulses, demonstrating potential future applications in quantum control and noisy intermediate-scale quantum computing.
Future planned enhancements include full integration into QuTiP's quantum control and QIP libraries, support for time-dependent bath parameters (e.g., timedependent chemical potential in fermionic systems), and support for solving problems with a system coupled to a fermionic and bosonic environment simultaneously within the HEOM description [74].For improved efficiency, further optimization of the fermionic solver construction (e.g., employment of the hermiticity relation discussed in [41] and [32]) is planned, as well as more robust and powerful parallel computing support, including GPU support.

Figure 6 .
Figure 6.(a) shows a fit of the spectral density J(ω) from Eq. (29) with a single over-damped Drude-Lorentz spectral density function, giving γ = 637 cm −1 and λ = 300 cm −1 .(b) shows a fit instead using two under-damped brownian motion spectral densities, giving α 2 1 = 0.11ω 3 c,1 , ωc,1 = 900 cm −1 , Γ1 = 800 cm −1 , and α 2 2 = 0.25ω 3 c,2 , ωc,2 = 500 cm −1 , Γ2 = 990 cm −1 .(c) shows the HEOM and RC results obtained for dynamics of several matrix elements on the two-site dimer model using the fit from (a), while (d) shows the same using the fit from (b).The RC and HEOM produce comparable results, and the choice of fitting function has little influence on the system dynamics in this case.(e) shows the entanglement between system and RC mode(s) as determined by the negativity using the single RC mode model for the fit in (a) and and two RC modes using the fit in (b).

Figure 8 .
Figure 8. Example of dynamical decoupling showing the coherence ρ01 of a two-level system coupled to an environment with a Drude-Lorentz spectral density.Here ∆ = 0, λ = 1 × 10 −4 Vf , γ = 1 × 10 −2 Vf , T = 0.1 Vf , where Vf is the amplitude of the fast drive.For the slow drive we choose V = 0.02 Vf , and in both cases we choose a pause period of ∆t Vf = 10.The lower figure, (b), shows the drive amplitudes as a function of time.In (a) we see that very slow pulses perform much worse than fast pulses, as expected due to the competition with the bath memory time.Note that here and in Fig.9we operate in a rotating frame with the drive, and assume that the system frequency is much larger than the drive amplitude ( V ) such that the rotating wave approximation is valid.The code for generating this figure can be found in example notebook 4[26].

Figure 9
Figure 9. Example of dynamical decoupling showing the coherence ρ01, after a total evolution time Tmax = 1000(1/ V ), of a two-level system coupled to an environment with a Drude-Lorentz spectral density versus γ, the width (inverse memory time) of the bath.Each run includes 100 dynamical decoupling pulses, and we compare equally spaced pulses (labelled DD in the figure legend) with the UDD scheme described in the text.The parameters here are chosen as ∆ = 0, λ = 0.001 V , 0.0005 V , 0.00005 V , and T = 0.05 V .There is a clear crossover for larger values of γ where the UDD scheme becomes less efficient than the standard equally spaced scheme (DD).The code for generating this figure can be found in example notebook 4[26].

Figure 10 .
Figure 10.A demonstration of the fitting of the spectral density of an Ohmic environment with exponential cut-off [Eq.(36) with s = 1] using kJ = 4 underdamped Brownian motion terms.We use bath parameters α = 3.25 and T = ωc.The resulting correlation functions are shown in (a) and (b), the original and fit spectral density is shown in (c), and the resulting power spectrum is shown in (d).In evaluating the correlation functions we use only one Matsubara term K = 1.The residual error seen in the real part is then compensated for using a terminator.The code for generating this figure can be found in example notebook 1d[26].

Figure 11 .
Figure 11.A demonstration of the fitting of the correlation functions of an Ohmic environment with exponential cut-off [Eq.(36) with s = 1] using kR = 3 and kI = 3 exponential terms terms.We use bath parameters α = 3.25 and T = ωc.The original and fit correlation functions are show in (a) and (b).The resulting spectral density and power spectrum are also shown in (c) and (d) respectively.The code for generating this figure can be found in example notebook 1d [26].

3 Figure 12 .
Figure 12.Dynamics of the excited state population ρ11(t) for a two-level system coupled to an Ohmic environment with exponential cut-off [Eq.(36) with s = 1] for α = 3.25, T = 0.5ωc, = 0, ∆ = 0.2ωc.We show results for fitting the spectral density with kJ = 1, and 3 underdamped Brownian motion terms, and for fitting the correlation functions with kr = kI = 1 and 3 terms.For the case where we fit the spectral density, we use N k = 1 and the standard terminator to approximate the other Matsubara terms.Here, and in Fig.13a HEOM truncation of Nc = 11 is used to get convergence.Examples of these fits used to generate the data in this figure are shown in Fig.10and Fig.11.The code for generating this figure can be found in example notebook 1d[26].
Figure 12.Dynamics of the excited state population ρ11(t) for a two-level system coupled to an Ohmic environment with exponential cut-off [Eq.(36) with s = 1] for α = 3.25, T = 0.5ωc, = 0, ∆ = 0.2ωc.We show results for fitting the spectral density with kJ = 1, and 3 underdamped Brownian motion terms, and for fitting the correlation functions with kr = kI = 1 and 3 terms.For the case where we fit the spectral density, we use N k = 1 and the standard terminator to approximate the other Matsubara terms.Here, and in Fig.13a HEOM truncation of Nc = 11 is used to get convergence.Examples of these fits used to generate the data in this figure are shown in Fig.10and Fig.11.The code for generating this figure can be found in example notebook 1d[26].

4 Figure 13
Figure 13.Difference in the dynamics as evaluated by the HEOM and that given by an analytical result for the coherence |ρ10(t) − ρ ana 10 | for a two-level system initially prepared with maximal coherence, and coupled to an Ohmic environment with exponential cut-off [Eq.(36) with s = 1] for α = 3.25, = 0, and ∆ = 0.The latter means that, unlike Fig.12, this result is integrable, and we can compare to an analytical result.(a) is for T = 0.5ωc, as in Fig.12, while (b) is for a lower temperature of T = 0.25ωc.In (a) we show results for fitting the spectral density with kJ = 1, and 3 underdamped Brownian motion terms, and for fitting the correlation functions with kr = kI = 1 and 3 terms.For the case where we fit the spectral density, we here use N k = 1 and the standard terminator to approximate the other Matsubara terms.In (b) the lower temperature means we have to both increase the number of spectral functions and employ more Matsubara terms to reach results comparable to the correlation function fitting results.

2 Figure 14 .
Figure 14.The analytical current and the result from the HEOM, using the Padé decomposition, for the single-impurity model, as a function of the bias voltage ∆µ, for η = 0.01 eV, T = 300 K, and W = 1 eV, and the impurity energy is = 0.3 eV.The code for generating this figure can be found in example notebook 5a [26].

Figure 15 .
Figure 15.Current versus the bias voltage for a single impurity coupled to a vibronic mode with = 0.3, eV Ω = 0.2 eV, λ = 0.12 eV, T = 300 K, W = 10 4 eV, and η = 0.01 eV.Note that increasing lmax tends to increase the current at large bias voltages, while increasing the other convergence parameter, the bosonic cut-off N , tends to decrease it, gradually giving a converged result in the vicinity of the red curve.The code for generating this figure can be found in example notebook 5b[26].

Table I .
A table comparing available HEOM implementations.The terminology used is as follows.Name: Abbreviated name of the implementation.Docs: "Yes" if there seemed to be sufficient documentation for a third-party to install and run the software on their own problems, without having to read the source code.Tests: "Yes" if the software included a suite of tests that cover the main features of the implementation and checks the resulting outputs.Time-dependent Hamiltonians: "Yes" if time-dependent Hamiltonians are supported.Baths: A description of which types of baths are supported.DL refers to Drude-Lorentz, UD refers to Under-damped Brownian motion, "Lor" refers to Lorentzian (for fermionic baths), and "fit" refers to fitting via generic decomposition of correlation functions.Compute: A description of what computational paradigms are supported.
n = (n R1 , n R2 , ..n RN R , n I1 , n I2 , ..., n IN I ) is a multi-index of integers n jk ∈ {0..N c },and N c is a cut-off parameter chosen for convergence.The state labelled by and we assume each independent bath is