Getting Started

This page walks you through a minimal TENSO simulation. After reading it you will understand the basic workflow and be ready to explore the examples and the code structure.

Conceptual overview

A TENSO simulation always follows the same four steps:

Define the bath correlation function $C(t)$ and decompose it into features
Define the system Hamiltonian $H_\textrm{S}$, coupling operator $Q_\textrm{S}$, and initial state $\rho_0$
Set propagation parameters (TTN topology, rank, depth, time grid)
Propagate the system and collect results

All four steps are handled by gen_bcf and system_multibath from tenso.prototypes.

The bath and its BCF decomposition

TENSO treats environments as a collection of harmonic oscillators coupled through the system-bath Hamiltonian $H_\textrm{SB} = Q_\textrm{S} \otimes X_\textrm{B}$. All bath information is encoded in the Bath Correlation Function (BCF):

$$C(t) = \langle \tilde{X}_\textrm{B}(t)\,\tilde{X}_\textrm{B}(0)\rangle$$

Using the residue theorem, the BCF is decomposed into $K$ complex decaying exponentials — the features:

$$C(t) = \sum_{k=1}^{K} c_k\, e^{\gamma_k t} \quad \text{and} \quad C^*(t) = \sum_{k=1}^{K} \bar{c}_k e^{\gamma_k t}$$

gen_bcf performs this decomposition from a spectral density $J(\omega)$ built from Drude–Lorentz (DL) and/or Brownian oscillator (BO) components:

$$J(\omega) = \sum_a J_\textrm{DL}^{(a)}(\omega) + \sum_b J_\textrm{BO}^{(b)}(\omega)$$

Drude–Lorentz — Ohmic bath with reorganization energy $\lambda$ and cutoff frequency $\omega_c$:

$$J_\textrm{DL}^{(a)}(\omega) = \frac{2\lambda}{\pi}\frac{\omega_c\,\omega}{\omega^2 + \omega_c^2}$$

Each DL component contributes one feature (a simple decaying exponential with timescale $\omega_c^{-1}$).

Brownian oscillator — discrete vibrational mode with reorganization energy $\lambda$, natural frequency $\omega_0$, and damping rate $\eta$:

$$J_\textrm{BO}^{(b)}(\omega) = \frac{4\lambda}{\pi}\frac{\eta\,\omega_0^2\,\omega}{(\omega^2 - \omega_0^2)^2 + 4\eta^2\omega^2}$$

Each BO component contributes two features (oscillatory-decaying exponentials). The total number of features $K$ also grows with the number of low-temperature correction (LTC) terms required to accurately capture the Bose–Einstein distribution $f_\textrm{BE}(\beta\omega)$ at finite temperature.

The total number of features $K$ counts as: 1 per DL component + 2 per BO component + 1 per LTC term. For the call below, $K = 1 + 2 + 1 = 4$.

from tenso.prototypes.bath import gen_bcf

bath = gen_bcf(
    re_d    = [540],    # λ for DL component (cm⁻¹)
    width_d = [70],     # ωc for DL component (cm⁻¹)
    freq_b  = [1663],   # ω0 for BO component (cm⁻¹)
    re_b    = [330],    # λ  for BO component (cm⁻¹)
    width_b = [4],      # η  for BO component (cm⁻¹)
    temperature          = 300,    # K
    decomposition_method = 'Pade', # 'Pade' or 'Matsubara'
    n_ltc                = 1,      # number of low-temperature correction terms
)

The system

The system is an $M$-level quantum system, where $M$ is the dimension of its Hilbert space. Each auxiliary density matrix (ADM) in the HEOM hierarchy has the same $M \times M$ dimension as $\rho_\textrm{S}(t)$. For the spin-boson model with energy gap $\varepsilon$ and tunneling $\Delta$:

$$H_S = \frac{\varepsilon}{2}\,\sigma_z + \Delta\,\sigma_x$$

The bath couples to the system through $Q_\textrm{S} = \sigma_x$, where the Pauli matrices are

$$\sigma_z = \begin{bmatrix}1 & 0 \\ 0 & -1\end{bmatrix}, \qquad \sigma_x = \begin{bmatrix}0 & 1 \\ 1 & 0\end{bmatrix}$$

Specify the system with three NumPy arrays:

sys_ham — the $M \times M$ Hamiltonian $H_\textrm{S}$ (complex128)
sys_op — the $M \times M$ coupling operator $Q_\textrm{S}$ (complex128)
init_rdo — the $M \times M$ initial reduced density operator $\rho_S(0)$

import numpy as np

eps   = 1500.0   # cm⁻¹  energy gap ε
delta =  300.0   # cm⁻¹  tunneling Δ

sys_ham = np.array([[ eps/2,  delta],
                    [ delta, -eps/2]], dtype=np.complex128)

# Q_S — system-bath coupling operator (σ_x)
sys_op  = np.array([[0.0, 1.0],
                    [1.0, 0.0]], dtype=np.complex128)

# Initial state: |↑⟩ (excited state)
wfn      = np.array([1.0, 0.0], dtype=np.complex128)
init_rdo = np.outer(wfn, wfn.conj())

Propagation parameters

Four numerical parameters control accuracy and cost:

Symbol	Meaning	Set by	Notes
$M$	System dimension	shape of `sys_ham`	Fixed by the physical model; $M=2$ for a qubit
$K$	BCF features	`gen_bcf` output	1 per DL + 2 per BO + 1 per LTC
$N$	HEOM depth	`dim`	Increase until converged
$R$	Bond rank	`rank`	Increase until converged

Without compression, storing the full EDO $|\Omega(t)\rangle$ requires $\mathcal{O}(M^2 N^K)$ memory — exponential in $K$. TENSO's TTN decomposition reduces this to $\mathcal{O}(M^2 R + KNR(N+R))$, which grows only polynomially with $K$.

The remaining propagation parameters are:

Parameter	Argument	Typical value	Effect
HEOM depth	`dim`	10–20	Sets $N_k$ for each of the $K$ bexciton ladders
Bond rank	`rank`	5–60+	Sets $R_s$ for each of the $K-1$ TTN bonds
TTN topology	`frame_method`	`'tree2'` or `'train'`	Balanced tree (`tree2`) minimizes average path from root to each bexciton to $\mathcal{O}(\log K)$; tensor train has path $\mathcal{O}(K)$
Propagation method	`ps_method`	`'ps1'`, `'ps2'` or `'vmf'`	PS2 adapts $R$ automatically; PS1 keeps $R$ fixed
Time step	`step_time`	0.05 fs	Integration step size $\Delta t$
End time	`end_time`	problem-dependent	Total propagation time in fs
Output file	`fname`	any string	Prefix for output files: `{fname}.dat.log`, `{fname}.debug.log`, `{fname}.pt`

frame_method='tree2' builds a balanced binary tree, minimizing the average path from the system root $A^{(0)}$ to each bexciton index to $\mathcal{O}(\log K)$, giving more compact TTNs for large $K$. frame_method='train' uses a tensor train (MPS) topology with path length $\mathcal{O}(K)$, which is simpler but less efficient for structured baths.

Running a simulation

system_multibath returns a propagator that you drive with a for loop:

from math import ceil
from tqdm import tqdm
from tenso.prototypes.heom import system_multibath

end_time = 100.0   # fs
dt       = 0.05    # fs

propagator = system_multibath(
    fname             = 'out',        # output prefix
    init_rdo          = init_rdo,
    sys_ham           = sys_ham,
    sys_ops           = [sys_op],
    bath_correlations = [bath],
    dim               = 20,           # HEOM truncation depth N_k
    end_time          = end_time,
    step_time         = dt,
    frame_method      = 'tree2',
    rank              = 20,
    stepwise_method   = 'simple',
    ps_method         = 'ps1',
)

progress_bar = tqdm(propagator, total=ceil(end_time / dt))
for _t in progress_bar:
    progress_bar.set_description(f'@{_t:.2f} fs')

When the loop finishes, three output files are written:

File	Contents
`out.dat.log`	Space-separated text (complex128). Columns: `t, ρ_00, ρ_01, ρ_10, ρ_11` (row-major flattening of $\rho_\textrm{S}$)
`out.debug.log`	Human-readable convergence diagnostics
`out.pt`	PyTorch checkpoint of the full TTN state at the final time step

Load and plot results:

import numpy as np
import matplotlib.pyplot as plt

data = np.loadtxt('out.dat.log', dtype=np.complex128)

t       = data[:, 0].real    # time (fs)
pop_exc = data[:, 1].real    # excited state population  [ρ_S(t)]_00  (|↑⟩ is index 0)
coh_eg  = data[:, 2]         # coherence ρ_eg = [ρ_S(t)]_01  (complex)
coh_ge  = data[:, 3]         # coherence ρ_ge = [ρ_S(t)]_10 = ρ_eg*  (complex)
pop_gnd = data[:, 4].real    # ground state population   [ρ_S(t)]_11

fig, axes = plt.subplots(1, 2, figsize=(8, 3))
axes[0].plot(t, pop_exc, label=r'$[\rho_\mathrm{S}]_{00}$')
axes[0].plot(t, pop_gnd, '--', label=r'$[\rho_\mathrm{S}]_{11}$')
axes[0].set(xlabel='Time (fs)', ylabel='Population')
axes[0].legend()
axes[1].plot(t, np.abs(coh_eg), color='tab:purple')
axes[1].set(xlabel='Time (fs)', ylabel=r'$|\rho_{eg}|$')
plt.tight_layout()
plt.show()

To reload the TTN state checkpoint:

import torch
state = torch.load('out.pt')

Checking convergence

TTN-HEOM converges to exact HEOM as $N$ and $R$ increase. Always verify convergence:

HEOM depth dim ($N_k$) — increase until $\rho_S(t)$ is stable
Bond rank rank ($R_s$) — increase until $\rho_S(t)$ is stable
Low-temperature corrections n_ltc — increase until long-time thermalization is correct

A typical convergence sweep:

from math import ceil
from tqdm import tqdm
from tenso.prototypes.heom import system_multibath
from tenso.prototypes.bath import gen_bcf

bath = gen_bcf(
    re_d=[540], width_d=[70],
    freq_b=[1663], re_b=[330], width_b=[4],
    temperature=300, decomposition_method='Pade', n_ltc=1,
)

end_time      = 100.0   # fs
dt            = 0.05    # fs
ranks         = [5, 10, 15, 20, 25, 32]
frame_methods = ['train', 'tree2']

for method in frame_methods:
    for rank in ranks:
        fname = f'{method}_rank{rank}'
        propagator = system_multibath(
            fname=fname, init_rdo=init_rdo,
            sys_ham=sys_ham, sys_ops=[sys_op],
            bath_correlations=[bath],
            dim=20, end_time=end_time, step_time=dt,
            frame_method=method, rank=rank,
            stepwise_method='simple', ps_method='ps1',
        )
        progress_bar = tqdm(propagator, total=ceil(end_time / dt))
        for _t in progress_bar:
            progress_bar.set_description(f'@{_t:.2f} fs')

Overlay results to verify convergence:

import numpy as np
import matplotlib.pyplot as plt

fig, ax = plt.subplots()

for method in ['train', 'tree2']:
    for rank in [5, 10, 20, 32]:
        data    = np.loadtxt(f'{method}_rank{rank}.dat.log', dtype=np.complex128)
        t       = data[:, 0].real
        pop_exc = data[:, 1].real
        ax.plot(t, pop_exc, label=f'{method} R={rank}')

ax.set(xlabel='Time (fs)', ylabel=r'$[\rho_\mathrm{S}]_{00}$')
ax.legend(fontsize=7)
plt.show()

Background Examples