Auxiliary Density Operators in HEOM

• Auxiliary Density Operators are mathematical constructs in the HEOM formalism that encode system-bath interactions and capture non-Markovian memory effects.
• They map environmental correlations into a hierarchical ladder of coupled evolution equations, enabling numerically exact simulations of complex quantum dynamics.
• Implementations like HierarchicalEOM.jl demonstrate their practical utility in modeling spectroscopy and dynamics in molecular aggregates and quantum dots.

Auxiliary Density Operators (ADOs) are a central mathematical construct in the hierarchical equations of motion (HEOM) formalism for simulating the nonequilibrium and dissipative dynamics of open quantum systems. By extending the dynamical description into a hierarchy of operators, ADOs enable a nonperturbative and numerically exact treatment of system-bath coupling, memory effects, and complex spectral densities, including those arising in systems with both bosonic and fermionic reservoirs. The auxiliary operators serve both as a compact representation of environmental correlations and as a computational tool by mapping the problem onto a ladder of coupled evolution equations, preserving information on bath-induced processes up to the truncation order.

1. Mathematical Foundation and Definition

In the HEOM approach, the combined system–environment Hamiltonian comprises the system Hamiltonian $\hat{H}_S$, a sum of bath Hamiltonians, and linear coupling terms $\sum_a \hat{V}_a \hat{X}_a$ where $\hat{X}_a$ represents collective bath coordinates. The reduced system dynamics are determined by tracing out environmental degrees of freedom, resulting in non-Markovian equations due to system-bath memory. Expansion of the bath correlation functions $C_a(t)=\langle \hat{X}_a(t)\hat{X}_a(0)\rangle$ as a sum of exponentials,

$C_a(t)=\sum_{k=0}^K c_{a,k}e^{-\gamma_{a,k}t},$

is a key step which allows for the introduction of ADOs indexed by non-negative integer vectors $n = \{n_{a,0}, n_{a,1}, \ldots, n_{a,K}\}$. The physical reduced density operator is $\rho_0(t)=\rho_S(t)$, with each higher ADO $\rho_n(t)$ encoding finer orders of system–bath interaction history (Seibt et al., 2018, Huang et al., 2023).

2. Hierarchy Structure and Recursion Relations

ADOs are arranged in a $K$-dimensional lattice, where tier $T:=|n|=\sum_{a,k} n_{a,k}$ groups all operators with the same total excitation number in the Matsubara decomposition. Each ADO $\rho_n$ evolves according to a coupled set of equations:

$$\partial_t \rho_n = -\left(i\mathcal{L}_S + \sum_{a,k} n_{a,k}\gamma_{a,k}\right)\rho_n - i\sum_{a,k}[V_a, \rho_{n+e_{a,k}}] - i\sum_{a,k} n_{a,k}\left[ c_{a,k}V_a\rho_{n-e_{a,k}} - \rho_{n-e_{a,k}}c_{a,k}^*V_a \right],$$

where $\mathcal{L}_S\cdot=[\hat{H}_S,\cdot]$ and $e_{a,k}$ is the unit vector in the $(a,k)$ direction. These couplings realize a structure analogous to creation and annihilation operators in bosonic mode occupation, with "hierarchy-up" and "hierarchy-down" mappings. Rescaling $$\tilde{\rho}_n = \rho_n/\sqrt{\prod_{a,k} n_{a,k}! |c_{a,k}|^{n_{a,k}}}$$ leads to a symmetric form where the coefficients directly mirror bath expansion amplitudes (Seibt et al., 2018).

3. Physical Interpretation and Analogy to Vibronic Basis

Each index $n_{a,k}$ in the ADO hierarchy quantifies the number of "quanta" associated with the $k$-th exponential of the $a$-th correlation function, mirroring vibrational mode excitations. The operators connecting different tiers act analogously to bosonic creation ($\sqrt{n_{a,k}+1}$ “up”) and annihilation ($\sqrt{n_{a,k}}$ “down”) operators. This correspondence connects ADOs to explicit vibrational bases in open quantum systems, where HEOM's bookkeeping over Matsubara terms effectively simulates vibronic progression without requiring an explicit enlarged basis (Seibt et al., 2018).

4. Extensions: Inclusion of Herzberg–Teller and Non-Condon Effects

Standard HEOM with ADOs is extended to model nonadiabatic Herzberg–Teller (HT) coupling and non-Condon effects by Taylor-expanding electronic couplings $J_{mn}(q)$ and transition dipoles $\mu_m(q)$ in vibrational coordinates $q_m$. This yields additional terms that couple different ADOs:

• First- and second-order HT coupling introduce new superoperators linking ADOs that differ by multiple indices, with prefactors determined by the HT derivatives and the bath expansion amplitudes.
• Non-Condon effects are incorporated by operator substitutions, introducing further cross-couplings in the ADO network. Such extensions increase the density of connections in the hierarchy, but preserve the general recursive structure. Consequently, coordinate-dependent electronic couplings and dipole moments can be captured entirely within the HEOM+ADO formalism, without explicit vibronic states (Seibt et al., 2018).

5. Hierarchy Truncation and Convergence Strategies

The multidimensional ADO structure is infinite, necessitating systematic truncation schemes:

• Tier cutoff: Limit total excitation number $|n|\leq N_\text{max}$.
• Importance filtering: Define an “importance” measure for each ADO, e.g.,

$$I(\rho^{(m,n,p)}_{j|q}) = \left|\prod_{r=1}^m \frac{\xi_{j_r}}{\mathrm{Re}[\chi_{j_r}]\,\sum_{x=1}^r \mathrm{Re}[\chi_{j_x}]}\right| \left|\prod_{w=1}^n \frac{\eta_{q_w}}{\mathrm{Re}[\gamma_{q_w}]\,\sum_{x=1}^w \mathrm{Re}[\gamma_{q_x}]}\right|.$$

ADOs with $I$ below a threshold $I_\text{th}$ are discarded. Empirically, $I_\text{th}\sim10^{-6}$--$10^{-7}$ is sufficient for convergence in typical scenarios. These schemes drastically reduce computational cost while maintaining numerical accuracy of system observables (Huang et al., 2023).

6. Computational Implementation and Practical Applications

The implementation of ADO-based HEOM, as exemplified in HierarchicalEOM.jl, involves:

• Encoding multi-index ADOs via a hierarchy dictionary associating flattened integer indices with multi-index data to maximize lookup and memory efficiency.
• Sparse assembly of the full HEOM Liouvillian ${\mathcal{M}}$ of dimension $d_\text{sys}^2 \times N_\text{ADO}$, where $N_\text{ADO}$ is the number of retained ADOs.
• Parallelized construction and solution of the coupled differential equations using Julia’s thread-based parallelism.
• Interfacing with high-performance ODE solvers and linear stationary-state routines to enable time propagation, stationary state search, and spectral computations.
• Direct extraction of physical observables—including the reduced density matrix, currents, and spectra—via targeted inspection or contraction of specific ADOs (e.g., first-tier fermionic ADOs provide current operators).

Table: Illustrative organization of the ADO hierarchy in HEOM

Tier	Multi-index example $n$	Operator encoded
0	$(0,0,\ldots,0)$	Physical reduced density operator $\rho_S(t)$
1	$(1,0,\ldots,0)$, $(0,1,0,\ldots)$	First-tier ADOs, encode leading bath memory
2	$(2,0,\ldots)$, $(1,1,0,\ldots)$, etc.	Second-tier ADOs, encode higher-order memory terms

HierarchicalEOM.jl provides a concrete realization of these ideas for both bosonic and fermionic hybrid baths (Huang et al., 2023).

7. Role in Optical Spectroscopy and Complex System Modeling

ADOs in HEOM have enabled detailed studies of linear absorption and nonlinear spectroscopy of vibronic dimers, demonstrating that bath-mediated effects, HT, and non-Condon corrections can be faithfully captured—even permitting calculation of effective temperature-dressed Huang–Rhys factors in the strong coupling limit. Calculated absorption spectra on the basis of ADOs show direct correspondence with results from explicit vibronic configurations, but with lower computational complexity due to the reduced operator formalism. The generalized ADO approach thus underpins accurate modeling of quantum dynamics in photosynthetic complexes, molecular aggregates, and quantum dots subjected to complex environments (Seibt et al., 2018).