TENSO Structure

TENSO is organized into four distinct layers of abstraction, each building on the one below it. The code structure spans these four layers plus a set of shared helper modules.

Overview

Layer 1 — Backend Interfaces

Module: tenso.backend

The tenso.backend module is the foundation upon which all computation in TENSO rests. It has three responsibilities: managing third-party dependencies, controlling global package settings (numerical precision, default device), and registering the integration methods available for time propagation. All higher layers interact with hardware and solvers exclusively through this module, which means that changing the ODE solver or hardware target requires touching only this one file.

Dependencies and design rationale

TENSO depends on three external libraries. The choice of each is deliberate:

PyTorch — primary tensor computation engine. PyTorch was chosen over NumPy-only or JAX-based alternatives for its hardware-agnostic tensor operations — the same code runs on CPU or GPU without modification, which is critical for the large tensor contractions that arise in TTN propagation. PyTorch's dynamic computation graph and automatic differentiation infrastructure also provide a flexible foundation for future extensions, such as gradient-based optimization of bath parameters. The backend exposes PyTorch's device management so that all TENSO objects (state tensors, operator matrices, propagator buffers) are consistently allocated on the same device.

NumPy — auxiliary numerics and interoperability. NumPy handles operations outside the main propagation loop — such as constructing spectral density grids, initializing sparse operator dictionaries, and interfacing with external analysis tools. It also serves as the interchange format when passing data to and from plotting libraries or other Python scientific packages. The backend ensures that NumPy arrays and PyTorch tensors can be converted without copying data where possible.

torchdiffeq — ODE integration. Rather than implementing ODE solvers from scratch, TENSO delegates time integration to torchdiffeq, a library of differentiable ODE solvers built on top of PyTorch. This gives access to high-quality adaptive-step solvers with error control, all operating natively on PyTorch tensors. The backend imports torchdiffeq and also implements any additional integration methods not covered by the library directly.

Global package settings

Beyond dependency management, tenso.backend enforces several global settings that affect all computations:

• Floating-point precision — the default dtype (torch.float64 for most quantum dynamics problems, where accumulated round-off in long propagations is a practical concern).
• Default device — CPU or GPU, set once and inherited by all subsequently constructed TENSO objects.
• Complex number support — TENSO works with complex-valued tensors throughout; the backend ensures that PyTorch's complex dtype handling is configured consistently.

Changing any of these settings after objects have been constructed is not supported; configure the backend before instantiating any TENSO objects.

Integration methods

Time propagation in TENSO reduces to integrating a large, structured ODE of the form $\dot{\Omega} = f(\Omega, t)$. The choice of solver has a significant effect on both accuracy and wall-clock time. The default is an adaptive Runge-Kutta solver, which works well for most problems. Stiff hierarchies — which arise when the bath relaxation timescales are much shorter than the system dynamics — may require implicit methods.

Layer 2 — TTN Decomposition

The general master equations addressed by TENSO take the form

where $|\Omega(t)\rangle$ is the Extended Density Operator (EDO) — a high-dimensional tensor that collects the physical reduced density matrix $\rho_\textrm{S}(t)$ together with the full hierarchy of auxiliary density matrices (ADMs). Concretely, the EDO is organized as

where $\vec{n} = (n_1, \ldots, n_K)$ indexes the bexcitonic occupation numbers and $\varrho_{\vec{0}}(t) = \rho_\textrm{S}(t)$ is the system's reduced density matrix. $\mathcal{L}(t)$ is the Liouvillian superoperator generating the dynamics. Layer 2 implements the efficient propagation of this equation without ever forming $|\Omega(t)\rangle$ or $\mathcal{L}(t)$ explicitly as dense objects.

TENSO addresses the $\mathcal{O}(M^2 N^K)$ exponential cost of the full EDO by decomposing it into a Tree Tensor Network (TTN) — a contraction of $K$ low-order core tensors $A^{(0)}(t), U^{(1)}(t), \ldots, U^{(K-1)}(t)$:

For bond dimensions $R_s = \mathcal{O}(R)$ and bexciton depth $N_k = \mathcal{O}(N)$, the space complexity of the TTN is $\mathcal{O}(M^2 R + KNR(N + R))$, which grows only polynomially with $K$. The root tensor $A^{(0)}(t)$ carries the system indices $i,j$ and the compressed index $a_1$, so the system's unitary dynamics is treated exactly while the bexciton influence is captured in compressed form. The semi-unitary core tensors $U^{(s)}(t)$ satisfy the orthogonality condition

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

The underlying graph of the TTN is a tree — there are no loops. This is the natural structure for both HEOM (where bexcitons are organized hierarchically) and ML-MCTDH. The absence of loops makes the key operations — tensor contractions, TDVP sweeps, singular value decompositions — computationally efficient. TENSO supports two canonical TTN topologies: the balanced tensor tree (minimizes the average path length between the system and each bexciton to $\mathcal{O}(\log K)$) and the tensor train / matrix product state form (linear chain, average path length $\mathcal{O}(K)$).

The Sum-of-Product generator

Even with a compact TTN representation of $|\Omega(t)\rangle$, applying $\mathcal{L}(t)$ is only efficient if $\mathcal{L}(t)$ itself has a structured form. In the bexcitonic formulation, the Liouvillian superoperator takes a natural Sum-of-Product (SoP) form:

where each product term consists of a component $h^>_m(t)$ acting on the ket index $|i\rangle$, a component $h^<_m(t)$ acting on the bra index $\langle j|$, and local operators $h^{(k)}_m$ acting on the $k$-th bexciton space $|n_k\rangle$. Each dissipator $\mathcal{D}_k$ in the bexcitonic master equation contributes five product terms, and the system Liouvillian $-iH_\textrm{S}^\times(t)$ contributes two more, giving $5K+2$ terms in total.

This factored structure means that applying one product term to the TTN requires only a sequence of small local contractions — one per node — rather than a global dense matrix-vector product. The SoP form is not an approximation: it follows exactly from the structure of the bexcitonic equations of motion.

The propagation pipeline

Layer 2 implements three objects that together form the propagation pipeline:

Frame  ──(topology)──►  Model  |Ω(t)⟩  =  Con(A⁽⁰⁾, U⁽¹⁾, …, U⁽ᴷ⁻¹⁾)
                                    │
                        SparseSPO   ℒ(t) = Σₘ h>ₘ ⊗ h<ₘ ⊗ h⁽¹⁾ₘ ⊗ … ⊗ h⁽ᴷ⁾ₘ
                                    │
                                    ▼
                      SparsePropagator ──► |Ω(t + Δt)⟩

A Frame defines the topology of the TTN (which nodes exist, how they connect, which bonds carry system vs. bexciton indices). A Model fills that topology with concrete tensor data — the root tensor $A^{(0)}(t)$ and the semi-unitary core tensors $\{U^{(s)}(t)\}$. A SparseSPO encodes the Liouvillian superoperator $\mathcal{L}(t)$ in SoP form. A SparsePropagator takes a Model and a SparseSPO and advances $|\Omega(t)\rangle$ in time.

tenso.state — State Representation

This module implements the data structure for $|\Omega(t)\rangle$, organized into two submodules: pureframe for the abstract graph infrastructure and puremodel for the concrete TTN vector.

tenso.state.pureframe — defines the graph structure of the tensor network independently of any numerical data. Separating topology from data allows the same Frame to be reused across multiple Model instances (e.g. during a basis change or when restarting from a checkpoint).

tenso.state.puremodel

tenso.operator.sparse — Generator & Propagator

This module encodes $\mathcal{L}(t)$ in SoP form and implements its action on a Model. Both classes work together: SparseSPO describes what operator to apply; SparsePropagator describes how and when to apply it.

What is the Dirac–Frenkel TDVP and why use it?

The Dirac–Frenkel Time-Dependent Variational Principle (TDVP) determines the optimal equations of motion for the core tensors $A^{(0)}(t)$ and $\{U^{(s)}(t)\}$ by requiring that the residual $\mathcal{L}(t)|\Omega(t)\rangle - \frac{d}{dt}|\Omega(t)\rangle$ is orthogonal to all allowed variations of $|\Omega(t)\rangle$ within the TTN ansatz:

This projection guarantees that the propagated state remains on the TTN manifold (bond dimensions stay bounded) while capturing the dynamics as accurately as possible for the given TTN ansatz. The equations of motion for the root tensor $A^{(0)}(t)$ and for each semi-unitary tensor $U^{(s)}(t)$ are derived from this principle together with the gauge condition that enforces semi-unitarity of the $U^{(s)}(t)$ throughout the propagation.

Propagation algorithms

Three propagation algorithms are implemented as methods of SparsePropagator. All three yield equivalent dynamics in the converged limit; numerical stability may vary with system parameters. The choice between them involves tradeoffs in memory, accuracy, and computational cost.

Direct Integration. All core tensor ODEs — for $A^{(0)}(t)$ and all $U^{(s)}(t)$ — are integrated simultaneously using the ODE solver configured in tenso.backend (default: adaptive RK45). This enables high-order time integration $\mathcal{O}(\Delta t^n)$ for an order-$n$ method and adaptive time stepping. Because the matrix $D^{(s)}(t)$ in the semi-unitary equation of motion is singular at $t = 0$ (for initially separable states), a regularization scheme is applied: singular values of $D^{(s)}$ below a threshold $\epsilon$ are replaced by $\epsilon$, introducing an error of only $\mathcal{O}(\epsilon^2)$. The singularity resolves after a brief initial transient (typically $<2$ fs).

PS1 — Fixed-rank Projector-Splitting. Based on a second-order symmetric Lie–Trotter splitting of $\mathcal{L}(t)$. Instead of propagating all tensors simultaneously, the algorithm 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. The forward and backward sweeps together constitute one time step $\Delta t$, achieving second-order $\mathcal{O}(\Delta t^3)$ Trotter accuracy. Bond dimensions remain constant throughout — PS1 is a fixed-rank method. It requires no regularization but is harder to parallelize than direct integration due to the sequential SVD structure.

PS2 — Adaptive-rank Projector-Splitting. Extends PS1 with a two-site move that combines adjacent core tensors into a fourth-order tensor, propagates it, and then decomposes it back via a truncated SVD with threshold $\epsilon'$. This allows bond dimensions to grow or shrink adaptively during propagation, with the rank at each bond controlled by the SVD truncation error. PS2 is the recommended starting point when the required ranks are not known in advance: it determines the appropriate computational resources automatically from the early-time dynamics.

Layer 3 — Master Equation Implementations

TENSO currently provides built-in implementations for two major quantum dynamics methods. Both are cast as special cases of the general TTN master equation and share an identical downstream interface: a Hierarchy object constructs the required SparsePropagator and Model for the given physical problem.

tenso.heom — Hierarchical Equations of Motion

HEOM is an exact, non-perturbative method for simulating the open quantum dynamics of driven quantum systems coupled to bosonic thermal baths. TENSO implements the bexcitonic generalization of HEOM, which provides a unified framework for many HEOM variants and admits arbitrary time-dependence in the system Hamiltonian $H_\textrm{S}(t)$.

Physical background. The influence of a thermal bath on the system is completely characterized by its Bath Correlation Function (BCF):

where $c_k$, $\bar{c}_k$, $\gamma_k \in \mathbb{C}$ are obtained by decomposing the spectral density $J(\omega)$ via Padé or Matsubara expansion of the Bose–Einstein factor $f_\textrm{BE}(\beta\omega)$. Each $k$ in this sum defines a feature of the bath. The number of features $K$ grows with the complexity of $J(\omega)$ and with the number of low-temperature correction terms required. The bexcitonic picture associates a fictitious bosonic quasiparticle — a bexciton of label $k$ — with each feature. The bexciton creation and annihilation operators $\hat{\alpha}^\dagger_k$, $\hat{\alpha}_k$ connect the ADMs in the hierarchy and give the bexcitonic master equation its ladder structure. The full bexcitonic equation of motion is

where $\mathcal{D}_k = \gamma_k\hat{\alpha}^\dagger_k\hat{\alpha}_k + (c_k Q_\textrm{S}^> - \bar{c}_k Q_\textrm{S}^<)\hat{z}_k^{-1}\hat{\alpha}^\dagger_k - Q_\textrm{S}^\times\hat{\alpha}_k\hat{z}_k$ is the dissipator associated with the $k$-th bath feature and $\hat{z}_k$ is an invertible metric operator. The standard HEOM is recovered as the specific case with number representation and $\hat{z}_k = i(\hat{\alpha}^\dagger_k\hat{\alpha}_k)^{-1/2}$. See the publications for the full derivation of the bexcitonic master equation.

Two submodules are provided depending on the number of baths:

Key classes (available in both eom and meom):

tenso.mctdh — MCTDH and ML-MCTDH

The Multiconfiguration Time-Dependent Hartree (MCTDH) method and its multilayer extension (ML-MCTDH) are implemented in tenso.mctdh. The three core design principles of TTN-HEOM — a Sum-of-Product generator, a TTN decomposition of the state, and the Dirac–Frenkel TDVP — are identical to those of ML-MCTDH. The key difference is that ML-MCTDH targets unitary dynamics of closed quantum systems (wavefunctions $|\psi(t)\rangle$ propagated by the Schrödinger equation), whereas TTN-HEOM targets non-unitary dissipative dynamics of open quantum systems (the EDO $|\Omega(t)\rangle$ propagated by the bexcitonic master equation). Both are implemented in TENSO through the same Model / SparsePropagator interface, with the physical differences encoded entirely in the SparseSPO and the initial conditions of the Model.

Layer 4 — Templates for Common Problems

Subpackage: tenso.prototypes

The tenso.prototypes subpackage provides ready-to-use interface templates for the most common physical problems. These functions abstract away all details of tensor network construction, hierarchy generation, and propagator setup.

Helper & Auxiliary Modules

Several auxiliary modules provide supporting capabilities shared across multiple layers.