Background

Before diving into the workflow, this page introduces the key concepts and variables used throughout TENSO's documentation. Understanding the Extended Density Operator and the exponential scaling problem it solves will help you make informed choices about TTN topology and propagation parameters.

The Extended Density Operator (EDO) and Bexcitons

HEOM simulates the open quantum dynamics of a quantum system coupled to a bosonic thermal bath by tracking not only the physical reduced density matrix $\rho_\textrm{S}(t)$ but also an entire hierarchy of auxiliary density matrices (ADMs) that encode the memory of the bath. In the bexcitonic formulation, these are organized into a single high-dimensional tensor — the Extended Density Operator (EDO):

$$|\Omega(t)\rangle = \sum_{\vec{n}}\, \varrho_{\vec{n}}(t)\,|\vec{n}\rangle$$

where $\vec{n} = (n_1, \ldots, n_K)$ is a multi-index of bexciton occupation numbers and $\varrho_{\vec{0}}(t) = \rho_\textrm{S}(t)$ is the physical reduced density matrix.

Each index $n_k$ runs from $0$ to $N_k - 1$, where $N_k$ is the truncation depth of the $k$-th bexciton ladder. The system itself lives in an $M$-dimensional Hilbert space, so each ADM $\varrho_{\vec{n}}(t)$ is an $M \times M$ matrix.

What are bexcitons?

Bexcitons are fictitious bosonic quasiparticles — one per bath feature — that capture the influence of the bath on the system. During the open quantum dynamics, bexcitons are created and destroyed until the system reaches thermal equilibrium. The bexcitonic creation and annihilation operators connect the different ADMs in the HEOM. The bexcitonic picture unifies many HEOM variants.

The full bexcitonic equation of motion for the EDO takes the form:

$$\frac{\mathrm{d}}{\mathrm{d}t}|\Omega(t)\rangle = \left(-iH_\textrm{S}^\times(t) + \sum_{k=1}^{K} \mathcal{D}_k\right)|\Omega(t)\rangle$$

where $\mathcal{D}_k$ is the dissipator associated with the $k$-th bath feature. The standard HEOM is recovered as a specific case. See the publications for the full derivation.

Variables at a glance

The four key numerical parameters in any TENSO simulation are:

Symbol Name Physical meaning How to set
$M$ System dimension Dimension of the system Hilbert space. Each ADM is an $M \times M$ matrix. $M = 2$ for a two-level system (qubit). Fixed by the physical model; determined by the shape of sys_ham.
$K$ Number of features Number of terms in the BCF decomposition $C(t) = \sum_{k=1}^{K} c_k e^{\gamma_k t}$. Each feature corresponds to one bexciton and one index of the EDO. Set by the bath model and temperature; output of gen_bcf. Equals 1 per DL + 2 per BO + 1 per LTC term.
$N$ HEOM depth (truncation) Maximum bexciton occupation $N_k - 1$ for each ladder. Controls how many levels of the hierarchy are retained. Must be increased until results are converged. Set via dim in the propagation call.
$R$ Bond rank (TTN compression) Dimension of the contracted bonds $\{a_s\}$ in the TTN. Controls compression quality. Must be increased until results are converged. Set via rank in the propagation call.
Both $N$ (dim) and $R$ (rank) are convergence parameters. Always increase them systematically until the observable $\rho_\textrm{S}(t)$ stops changing. See Checking convergence for a practical sweep strategy.

Why tensor networks?

The full EDO tensor has elements $[\Omega(t)]_{i,j,n_1\cdots n_K}$, where $i,j$ index the $M$-dimensional system Hilbert space and each $n_k$ indexes a bexciton ladder of depth $N_k$. Storing the EDO without compression requires

$$\mathcal{O}(M^2 N^K)$$

memory — exponential in $K$. For the thymine spectral density from the TENSO paper ($K = 20$, $N = 20$, $M = 2$), the uncompressed EDO would require roughly 6.7 ronnabytes — far beyond any present computer.

The exponential wall

Even with moderate parameters — $K = 10$ features, $N = 10$ truncation depth, $M = 2$ system states — the uncompressed EDO requires $2^2 \times 10^{10} \approx 4 \times 10^{10}$ complex numbers, or about 640 GB. TTN compression makes this tractable.

TENSO's Tree Tensor Network (TTN) decomposition replaces the full EDO with a product of $K$ low-order core tensors $A^{(0)}(t), U^{(1)}(t), \ldots, U^{(K-1)}(t)$, reducing the storage to

$$\mathcal{O}(M^2 R + KNR(N + R))$$

which grows only polynomially with $K$. The first term $M^2 R$ is the cost of the root tensor $A^{(0)}(t)$ that carries the uncompressed system indices $i, j$; the second term $KNR(N + R)$ is the cost of the $K$ semi-unitary core tensors $U^{(s)}(t)$.

The TTN decomposition of the EDO is written explicitly as:

$$[\Omega(t)]_{ij n_1\cdots n_K} = \sum_{a_1 \cdots a_{K-1}} A^{(0)}_{ij a_1}(t)\; U^{(1)}_{a_1\beta_1\gamma_1}(t)\;\cdots\; U^{(K-1)}_{a_{K-1}\beta_{K-1}\gamma_{K-1}}(t)$$

The semi-unitary core tensors $U^{(s)}(t)$ satisfy the orthogonality condition

$$\sum_{\beta\gamma} [U^{(s)}(t)]^*_{a'_s,\beta\gamma}\,[U^{(s)}(t)]_{a_s,\beta\gamma} = \delta_{a'_s a_s}$$

which is enforced throughout the propagation and is essential for numerical stability.

For a deeper look at the mathematical structure of the TTN decomposition and how the Liouvillian superoperator is represented in Sum-of-Product form, see the TENSO Structure page.

TTN topologies

TENSO supports two canonical topologies for the TTN, selected via the frame_method parameter:

Balanced tensor tree

frame_method='tree2' — Organizes the $K$ bexciton indices into a balanced binary tree rooted at $A^{(0)}$. The average path length from the system root to any bexciton index is $\mathcal{O}(\log K)$, which minimizes the compression cost for large $K$. Recommended for most structured baths.

Tensor train (MPS)

frame_method='train' — Arranges the bexciton indices in a linear chain. The average path length is $\mathcal{O}(K)$, which is simpler but typically less efficient for structured baths with many features. Easier to interpret and can be faster for very small $K$.

For most structured baths, tree2 is the recommended starting point. Both topologies converge to the same exact dynamics as $R \to \infty$. When comparing methods, always run both to verify consistency. Although TENSO currently has these two implemented structures, the user can define its own custom topologies always ensuring that the TN is a tree.
Property tree2 (balanced tree) train (tensor train)
Average path length $\mathcal{O}(\log K)$ $\mathcal{O}(K)$
Best for Large $K$, structured baths Small $K$, simple baths
Recommended Yes, default choice For comparison and small systems

Propagation strategies

TENSO implements three propagation algorithms, selected via the ps_method parameter. They can be switched on-the-fly during a simulation:

Direct Integration — ps_method='vmf'

All core tensors are integrated simultaneously using an adaptive ODE solver (default: RK45, i.e., dopri5). Provides fixed rank and high-order time accuracy $\mathcal{O}(\Delta t^n)$ for an order-$n$ solver. Requires regularization near $t = 0$ for initially separable states, but the singularity resolves after a brief initial transient period (typically <2 fs).

PS1 — Fixed-rank Projector-Splitting — ps_method='ps1'

Based on a second-order symmetric Lie–Trotter splitting of the Liouvillian. Sweeps through the TTN sequentially: at each split-step, the root is moved to the next tensor via SVD and that tensor is propagated in isolation. Forward and backward sweeps constitute one time step $\Delta t$, achieving $\mathcal{O}(\Delta t^3)$ Trotter accuracy. Bond dimensions remain constant throughout — no regularization needed.

PS2 — Adaptive-rank Projector-Splitting — ps_method='ps2'

Extends PS1 with a two-site move that combines adjacent core tensors, propagates the combined tensor, and decomposes it back via truncated SVD with threshold $\epsilon'$. This allows bond dimensions to grow or shrink adaptively. Recommended when the required ranks are not known in advance: PS2 determines the appropriate computational resources automatically from the early-time dynamics.

Property Direct Integration (vmf) PS1 (ps1) PS2 (ps2)
Rank during propagation Fixed Fixed Adaptive
Time accuracy $\mathcal{O}(\Delta t^n)$, solver-dependent $\mathcal{O}(\Delta t^3)$ (Trotter) $\mathcal{O}(\Delta t^3)$ (Trotter)
Regularization needed Yes (near $t = 0$) No No
Parallelizable Yes (all tensors simultaneously) Difficult (sequential SVD) Difficult (sequential SVD)
Best for Bulk dynamics; high-order time steps Known ranks; robust propagation Unknown ranks; automatic compression
A practical mixed strategy: begin with ps2 (which avoids the initial singularity and determines appropriate ranks automatically), then switch to direct integration (vmf) once adaptive ranks have stabilized. This achieves converged results with lower total computational cost.
For full mathematical details of the TDVP equations of motion, the Dirac–Frenkel variational principle, and algorithm derivations, see the TENSO Structure page and the publications.
Getting Started TENSO Structure Examples