Finite-Temperature Tensor Networks

• Finite-temperature tensor network methods are computational frameworks that represent and evolve thermal states via compressed tensor ansätze for efficient evaluation of thermal observables.
• They encompass diverse approaches such as direct operator evolution, purification, and vectorization, each balancing accuracy, computational cost, and representational constraints.
• Key challenges include managing Trotter and truncation errors, ensuring positivity, and performing accurate environment contractions, with tailored strategies for each method.

Finite-temperature tensor network methods are numerical schemes that represent thermal objects such as the Gibbs operator $\rho(\beta)=e^{-\beta H}/Z$, its purification, its vectorization, or an equivalent transfer network by compressed tensor-network ansätze, so that thermal observables $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$ can be evaluated in regimes where exact diagonalization is restricted and, in many relevant cases, quantum Monte Carlo is either unavailable or used mainly as a benchmark (Motoyama et al., 14 Jan 2025, Zhang et al., 2024). In current usage, the term covers direct tensor-product-operator evolution on infinite lattices, purification-based PEPS and iPEPS, vectorized PEPS with stochastic reconfiguration, thermal tensor renormalization of partition-function networks, chain-mapping methods for finite-temperature open quantum systems, and, in a distinct machine-learning sense, matrix-product-state classifiers endowed with a temperature layer (Kshetrimayum et al., 2018, Qian et al., 2024, Ueda et al., 7 Aug 2025, Lacroix et al., 2024, Lin et al., 2021).

1. Formal objects and tensor-network representations

For a time-independent Hamiltonian $\mathcal H$, the thermal state is the density operator

$$\rho(\beta)=\frac{e^{-\beta\mathcal H}}{Z(\beta)},\qquad Z(\beta)=\mathrm{Tr}\,e^{-\beta\mathcal H},$$

with $\beta=1/T$ (Motoyama et al., 14 Jan 2025). Finite-temperature tensor network methods differ primarily in which object is represented directly. One class represents $\rho(\beta)$ itself as a tensor product operator or projected entangled pair operator with two physical indices per site; another represents a purification $|\psi_\beta\rangle$ in an enlarged Hilbert space satisfying $$\rho_\beta=\mathrm{Tr}_Q|\psi_\beta\rangle\langle\psi_\beta|$$; a third vectorizes $\rho$ into $|\rho\rangle_\sharp$ and evolves it under a doubled-space generator such as $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$0 (Motoyama et al., 14 Jan 2025, Qian et al., 2024, Zhang et al., 2024). Thermal partition functions can also be treated directly as Euclidean tensor networks and coarse-grained by TRG, HOTRG, or TNR-type algorithms, rather than by explicit density-operator evolution (Kuramashi et al., 2018, Ueda et al., 7 Aug 2025).

This suggests a useful taxonomy of the field.

Representation	Core construction	Representative papers
Direct operator ansatz	iTPO/PEPO for $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$1	(Motoyama et al., 14 Jan 2025, Kshetrimayum et al., 2018)
Purification	$$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$2	(Qian et al., 2024, Sinha et al., 2022, Czarnik et al., 2019, Kadow et al., 2023)
Vectorization	$$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$3 in doubled space	(Zhang et al., 2024, Han et al., 24 Mar 2025)
Transfer / renormalization network	Partition-function or PEPO column coarse-graining	(Kuramashi et al., 2018, Kuramashi et al., 2018, Ueda et al., 7 Aug 2025, Wang et al., 2022)
Open-system embedding	Finite-$$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$4 bath encoded by chain mapping or mesoscopic reservoirs	(Lacroix et al., 2024, Brenes et al., 2019)

In direct operator methods, the local tensor carries bra and ket physical indices and several virtual legs. TeNeS-v2, for example, represents $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$5 as an infinite tensor product operator on an infinite square lattice, with local rank-6 tensor $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$6 and bond dimension $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$7 (Motoyama et al., 14 Jan 2025). In purification-based approaches, each site carries both a physical and an ancilla degree of freedom, and the infinite-temperature state is a product of maximally entangled local physical-ancilla pairs (Qian et al., 2024, Sinha et al., 2022, Kadow et al., 2023). In vectorized approaches, the operator is reinterpreted as a state on a doubled local space of dimension $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$8, so that imaginary-time evolution becomes a ground-state-like problem in enlarged Hilbert space (Zhang et al., 2024, Han et al., 24 Mar 2025).

The same phrase is also used outside many-body thermodynamics. In supervised learning, the “finite temperature tensor network” is an MPS classifier with a temperature layer $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$9, where the loss is interpreted as a thermal average rather than a simple sum over samples (Lin et al., 2021). That usage is formally related by the Boltzmann factor but is not a Gibbs-state solver.

2. Principal algorithmic constructions

A common construction is cooling from infinite temperature by imaginary-time evolution. In direct operator evolution, one starts from $\mathcal H$0 and applies Trotterized gates to bra and ket indices, truncating back to bond dimension $\mathcal H$1 after each step. In TeNeS-v2, one update with step size $\mathcal H$2 advances the physical inverse temperature by $\mathcal H$3,

$\mathcal H$4

with local gates $\mathcal H$5 and a simple-update truncation on the iTPO (Motoyama et al., 14 Jan 2025). The earlier annealing algorithm for two-dimensional thermal states follows the same logic in a PEPO language: begin from the infinite-temperature identity, apply imaginary-time gates in many small $\mathcal H$6 steps, and compress after each gate by simple update (Kshetrimayum et al., 2018).

Purification methods evolve a pure state in an enlarged Hilbert space. The standard form is

$\mathcal H$7

where $\mathcal H$8 purifies the infinite-temperature identity (Qian et al., 2024, Czarnik et al., 2019). In iPEPS implementations for the Kitaev-Heisenberg model and the infinite-square-lattice Hubbard model, the imaginary-time propagator is Trotterized into nearest-neighbor gates, each gate temporarily enlarges a bond, and a local variational truncation restores the target bond dimension (Czarnik et al., 2019, Sinha et al., 2022). The Hubbard work uses a fermionic neighborhood tensor update, while the Kitaev-Heisenberg work uses exact environment full update; both operate directly in the thermodynamic limit (Sinha et al., 2022, Czarnik et al., 2019). Isometric tensor network states provide a related purification framework in which an orthogonality column and row play the role of a two-dimensional canonical gauge; imaginary-time evolution is then implemented by TEBD$\mathcal H$9 with Moses Move splitting of double columns (Kadow et al., 2023).

Vectorization-based methods replace $$\rho(\beta)=\frac{e^{-\beta\mathcal H}}{Z(\beta)},\qquad Z(\beta)=\mathrm{Tr}\,e^{-\beta\mathcal H},$$0 by $$\rho(\beta)=\frac{e^{-\beta\mathcal H}}{Z(\beta)},\qquad Z(\beta)=\mathrm{Tr}\,e^{-\beta\mathcal H},$$1 and treat thermal evolution as a variational optimization in doubled space. In the scalable PEPS approach of Zhang et al., the vectorized thermal state satisfies

$$\rho(\beta)=\frac{e^{-\beta\mathcal H}}{Z(\beta)},\qquad Z(\beta)=\mathrm{Tr}\,e^{-\beta\mathcal H},$$2

and the parameters are updated by stochastic reconfiguration rather than by local Trotter-gate compression (Zhang et al., 2024). The Rydberg-array study adopts the same vectorization principle and uses stochastic reconfiguration on finite PEPS to access thermal criticality in long-range interacting systems (Han et al., 24 Mar 2025).

A different algorithmic branch works directly with partition-function networks. The three-dimensional $$\rho(\beta)=\frac{e^{-\beta\mathcal H}}{Z(\beta)},\qquad Z(\beta)=\mathrm{Tr}\,e^{-\beta\mathcal H},$$3 gauge-theory study rewrites the Euclidean partition function as a local tensor network with constraint tensors $$\rho(\beta)=\frac{e^{-\beta\mathcal H}}{Z(\beta)},\qquad Z(\beta)=\mathrm{Tr}\,e^{-\beta\mathcal H},$$4 and plaquette tensors $$\rho(\beta)=\frac{e^{-\beta\mathcal H}}{Z(\beta)},\qquad Z(\beta)=\mathrm{Tr}\,e^{-\beta\mathcal H},$$5, then applies HOTRG along the temporal direction and TRG in space (Kuramashi et al., 2018). The recent TTNR construction instead builds a PEPO for $$\rho(\beta)=\frac{e^{-\beta\mathcal H}}{Z(\beta)},\qquad Z(\beta)=\mathrm{Tr}\,e^{-\beta\mathcal H},$$6 by a cluster expansion, linearly stacks these layers into a column PEPO, traces physical indices to obtain a two-dimensional classical tensor network, and then coarse-grains that network with a global TNR scheme (Ueda et al., 7 Aug 2025). TN tailoring takes yet another route: start from the tensor network representing the zero-temperature partition function, “scissor” a finite segment from the infinite network, “stitch” it to impose periodic boundary conditions in imaginary time, and variationally fine-tune the boundary MPS to optimize the finite-$$\rho(\beta)=\frac{e^{-\beta\mathcal H}}{Z(\beta)},\qquad Z(\beta)=\mathrm{Tr}\,e^{-\beta\mathcal H},$$7 free energy (Wang et al., 2022).

For open quantum systems, finite temperature is often encoded in the environment rather than in a Gibbs-state tensor network. T-TEDOPA replaces a true finite-$$\rho(\beta)=\frac{e^{-\beta\mathcal H}}{Z(\beta)},\qquad Z(\beta)=\mathrm{Tr}\,e^{-\beta\mathcal H},$$8 bosonic bath by an equivalent zero-temperature chain with a temperature-dependent spectral density $$\rho(\beta)=\frac{e^{-\beta\mathcal H}}{Z(\beta)},\qquad Z(\beta)=\mathrm{Tr}\,e^{-\beta\mathcal H},$$9, after which the combined system and chain are evolved as an MPS or TTN by TDVP (Lacroix et al., 2024). The mesoscopic-reservoir method for thermal machines instead imposes Fermi-Dirac occupations through Lindblad damping of lead modes and evolves the resulting superfermion MPS to the non-equilibrium steady state (Brenes et al., 2019).

3. Contraction, environments, and computational structure

The central numerical bottleneck is contraction of the thermal network or, equivalently, construction of an effective environment for local updates and observables. Direct iTPO methods have a notable advantage: finite-temperature expectation values are single-layer contractions. In TeNeS-v2, $\beta=1/T$0 or $\beta=1/T$1 is reduced by contracting the bra and ket indices locally, producing a single-layer rank-4 tensor network that is then approximated by CTMRG. The paper emphasizes that finite-$\beta=1/T$2 iTPO contractions are single-layer, whereas ground-state iTPS expectations require a double-layer network, making finite-$\beta=1/T$3 calculations substantially cheaper at the same $\beta=1/T$4 (Motoyama et al., 14 Jan 2025). For $\beta=1/T$5, the reported CTMRG cost is $\beta=1/T$6 and the memory is $\beta=1/T$7 (Motoyama et al., 14 Jan 2025).

Purification and vectorization usually recover a doubled local space and therefore a more expensive contraction problem. Exact-environment full update for the honeycomb Kitaev-Heisenberg model avoids part of this cost by dynamically mapping the hexagonal lattice to a rhombic lattice of quadruple tensors, so that CTMRG only sees uniform $\beta=1/T$8 bonds even though the updated physical bond temporarily grows to $\beta=1/T$9 (Czarnik et al., 2019). The Hubbard purification algorithm uses neighborhood tensor update: a finite cluster around the updated bond is contracted exactly to build a Hermitian, non-negative metric tensor $\rho(\beta)$0, and the reduced matrices $\rho(\beta)$1, $\rho(\beta)$2 are optimized by alternating linear solves in that metric (Sinha et al., 2022).

Vectorized PEPS with stochastic reconfiguration shifts the computational burden from local environment optimization to Monte Carlo estimation of a metric matrix $\rho(\beta)$3 and a force vector $\rho(\beta)$4 in parameter space,

$\rho(\beta)$5

with matrix-vector products evaluated without explicitly forming the full $\rho(\beta)$6 (Zhang et al., 2024). In the Rydberg implementation, this makes long-range $\rho(\beta)$7 interactions natural, because the update is variational and global rather than gate-local (Han et al., 24 Mar 2025).

Thermal renormalization methods recast the environment problem as coarse-graining of a classical network. The global TNR scheme defines a local tensor difference $\rho(\beta)$8 inside a $\rho(\beta)$9 block and minimizes

$|\psi_\beta\rangle$0

where $|\psi_\beta\rangle$1 is obtained from CTMRG and $|\psi_\beta\rangle$2 in the reported simulations (Ueda et al., 7 Aug 2025). This augments the usual local truncation norm by an explicit environment-weighted term and is used inside TTNR after the finite-temperature PEPO column has been compressed to a two-dimensional network (Ueda et al., 7 Aug 2025).

Isometric tensor networks supply a different contraction paradigm. Because tensors away from the orthogonality hypersurface satisfy isometric constraints, contraction collapses exactly onto the orthogonality row and column, giving a two-dimensional analogue of canonical MPS contraction (Kadow et al., 2023). TN tailoring uses yet another environment representation: the dominant left and right boundary MPS of the zero-temperature transfer MPO define a variational boundary for the finite-height tailored network, and the free energy per site is optimized directly as

$|\psi_\beta\rangle$3

with respect to the boundary tensors (Wang et al., 2022).

4. Numerical errors, physical constraints, and recurring controversies

Across the literature, three approximations recur: Trotter or cluster-expansion error in imaginary time, truncation error from finite bond dimension, and environment error from approximate contraction. TeNeS-v2 states these explicitly for the iTPO algorithm and adds a fourth issue specific to direct operator methods: the approximate iTPO is not guaranteed to be positive semidefinite, so one may observe unphysical energies below the true ground-state energy; increasing the CTM bond dimension $|\psi_\beta\rangle$4 alone does not repair this, whereas increasing the iTPO bond dimension $|\psi_\beta\rangle$5 can (Motoyama et al., 14 Jan 2025). This is one of the main controversies in direct-operator thermal PEPS: single-layer contractions are cheaper, but positivity is not built in.

Purification addresses that particular issue because the density matrix is obtained by tracing an explicit pure state in an enlarged Hilbert space, but it pays with a doubled local description and generally double-layer contractions (Motoyama et al., 14 Jan 2025, Qian et al., 2024). Clifford-augmented TDVP modifies this tradeoff by transferring stabilizer-compatible entanglement to a Clifford circuit, so that only non-stabilizerness remains in the MPS. In the benchmarks reported, CA-TDVP with $|\psi_\beta\rangle$6 reaches approximately the same energy accuracy as ordinary TDVP with $|\psi_\beta\rangle$7 on the $|\psi_\beta\rangle$8 $|\psi_\beta\rangle$9–$$\rho_\beta=\mathrm{Tr}_Q|\psi_\beta\rangle\langle\psi_\beta|$$0 Heisenberg model, corresponding to an $$\rho_\beta=\mathrm{Tr}_Q|\psi_\beta\rangle\langle\psi_\beta|$$1 cost reduction that the paper estimates as roughly a factor of $$\rho_\beta=\mathrm{Tr}_Q|\psi_\beta\rangle\langle\psi_\beta|$$2 (Qian et al., 2024).

At criticality, finite bond dimension and finite environment dimension induce effective finite correlation lengths. The annealing PEPO study notes that finite $$\rho_\beta=\mathrm{Tr}_Q|\psi_\beta\rangle\langle\psi_\beta|$$3 can mimic finite-size behavior in superfluid regimes and must be scaled carefully near critical points (Kshetrimayum et al., 2018). The isoTNS study makes the limitation more explicit: finite-$$\rho_\beta=\mathrm{Tr}_Q|\psi_\beta\rangle\langle\psi_\beta|$$4 isoTNS support exponentially decaying correlations, so long-range critical correlations cannot be represented exactly; energy and correlation errors peak near the thermal critical point of the two-dimensional Ising model even when the exact non-isometric tensor network is known (Kadow et al., 2023). The Rydberg study adds a separate caveat: for noncommuting order parameters, the susceptibility entering finite-size scaling was approximated by the classical fluctuation-dissipation form rather than the full Kubo formula, because the latter was too expensive in their PEPS framework (Han et al., 24 Mar 2025).

Open-system chain-mapping methods introduce their own error hierarchy. In T-TEDOPA, accuracy is controlled by chain length $$\rho_\beta=\mathrm{Tr}_Q|\psi_\beta\rangle\langle\psi_\beta|$$5, bond dimension $$\rho_\beta=\mathrm{Tr}_Q|\psi_\beta\rangle\langle\psi_\beta|$$6, local bosonic cutoff $$\rho_\beta=\mathrm{Tr}_Q|\psi_\beta\rangle\langle\psi_\beta|$$7, and time step $$\rho_\beta=\mathrm{Tr}_Q|\psi_\beta\rangle\langle\psi_\beta|$$8; at finite temperature, the rule of thumb quoted for avoiding chain-end reflections is $$\rho_\beta=\mathrm{Tr}_Q|\psi_\beta\rangle\langle\psi_\beta|$$9, rather than $\rho$0 at zero temperature (Lacroix et al., 2024). The Bethe-lattice thermal TTN shows that even at a continuous thermal phase transition the correlation length can remain finite, reaching the lattice-specific maximum $\rho$1; this is a geometry-dependent result rather than a generic feature of finite-temperature tensor networks (Qu et al., 2019).

5. Benchmark problems and physical results

Two-dimensional Ising-type models remain the standard calibration ground. The annealing PEPO algorithm reproduces the finite-temperature phase transition of the Ising model on an infinite square lattice in remarkable agreement with the exact solution and then extends to hard-core and soft-core Bose-Hubbard models, where it captures temperature-driven suppression of the superfluid regime (Kshetrimayum et al., 2018). TeNeS-v2 benchmarks the transverse-field Ising model at $\rho$2 with $\rho$3 and $\rho$4: the $\rho$5 energy density deviates noticeably from QMC, while $\rho$6 and $\rho$7 almost coincide with QMC; at $\rho$8, energy density, specific heat, and $\rho$9 are all reported in good agreement with QMC (Motoyama et al., 14 Jan 2025).

Finite-temperature PEPS has also been used to extract universality data. For square-lattice Rydberg arrays at $|\rho\rangle_\sharp$0 and $|\rho\rangle_\sharp$1, the checkerboard transition yields $|\rho\rangle_\sharp$2 and finite-size-scaling estimates $|\rho\rangle_\sharp$3, $|\rho\rangle_\sharp$4, and $|\rho\rangle_\sharp$5, consistent with the two-dimensional Ising universality class; by contrast, the striated phase produces effective exponents $|\rho\rangle_\sharp$6 and $|\rho\rangle_\sharp$7 that the authors do not assign to a known universality class and attribute possibly to finite-size effects (Han et al., 24 Mar 2025). TTNR reaches the same style of conclusion from a transfer-matrix viewpoint: for the two-dimensional transverse-field Ising model, the thermal critical line yields central charge $|\rho\rangle_\sharp$8, and at one sample point $|\rho\rangle_\sharp$9 the reported value is $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$00 after 12 RG steps (Ueda et al., 7 Aug 2025).

Strongly correlated fermions illustrate the reach of thermodynamic-limit finite-temperature iPEPS. The infinite-square-lattice Hubbard study uses U(1)$$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$01U(1) symmetry and bond dimensions up to $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$02, reaching temperatures down to $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$03 times the hopping rate. It reports spin and charge correlators free of boundary effects, evidence that mobile holes disrupt the antiferromagnetic background in a slightly doped system, signatures of hole-doublon pairs near half filling, hole-hole repulsion upon doping, and specific heat in the slightly doped regime (Sinha et al., 2022). On finite clusters, the scalable vectorized PEPS method benchmarks the half-filled $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$04 Hubbard model against DQMC and reproduces the two-peak specific heat structure, double occupancy, and next-nearest-neighbor spin correlations (Zhang et al., 2024).

Frustrated magnets and spin liquids supply sterner tests. The finite-temperature iPEPS study of the Kitaev-Heisenberg model estimates $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$05 in the stripy regime and $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$06 in the antiferromagnetic regime, with duality mapping these to ferromagnetic and zigzag counterparts (Czarnik et al., 2019). In the pure Kitaev limit, the same method accurately recovers the high-temperature crossover to spin ordering and qualitatively the lower-temperature crossover to flux ordering, benchmarked against Monte Carlo (Czarnik et al., 2019).

Finite-temperature tensor networks also appear outside conventional two-dimensional equilibrium lattice thermodynamics. The three-dimensional $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$07 gauge-theory study uses HOTRG and TRG on the Euclidean partition function, reaches spatial sizes up to $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$08 for $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$09, and extracts $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$10 values $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$11, $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$12, and $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$13, consistent with the Svetitsky–Yaffe conjecture and two-dimensional Ising universality (Kuramashi et al., 2018). The Bethe-lattice XXZ work obtains a finite-temperature phase diagram with ferromagnetic, XY, antiferromagnetic, and paramagnetic phases, continuous thermal transitions at $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$14, first-order lines at $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$15, and Goldstone-like transverse modes characterized by the maximal correlation length $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$16 (Qu et al., 2019). For open quantum systems, MPSDynamics.jl benchmarks pure dephasing and driven spin-boson dynamics at zero and finite temperature, while the mesoscopic-reservoir method for thermal machines reproduces Landauer–Büttiker power and efficiency in the quadratic limit and then accesses interaction-enhanced power in a three-site heat engine beyond that limit (Lacroix et al., 2024, Brenes et al., 2019).

6. Scope, reinterpretations, and current directions

Finite-temperature tensor network method is therefore not a single algorithm but a methodological family whose members differ in the represented object, the contraction strategy, and the physical question being asked. Direct iTPO and PEPO schemes favor single-layer contractions and practical access to comparatively large bond dimensions in two dimensions, but they may violate positivity at finite $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$17 (Motoyama et al., 14 Jan 2025). Purification retains positivity by construction and fits naturally with TEBD or TDVP, but it usually incurs doubled local spaces and heavier contractions (Qian et al., 2024, Sinha et al., 2022). Vectorization plus stochastic reconfiguration treats finite-temperature evolution as a global variational problem and is particularly natural for long-range interactions and frustrated models (Zhang et al., 2024, Han et al., 24 Mar 2025). Thermal renormalization methods such as TTNR shift emphasis from order parameters and critical exponents toward direct extraction of CFT data, proposing central charge and transfer-matrix spectra as an alternative route to locating thermal phase transitions (Ueda et al., 7 Aug 2025).

Several recent proposals aim to bypass the usual low-temperature bottlenecks. Clifford augmentation reduces the effective entanglement burden of purification-based TDVP by moving stabilizer-compatible entanglement into a cheaply simulable Clifford circuit (Qian et al., 2024). TN tailoring reverses the standard annealing logic by starting from the zero-temperature tensor network and extracting finite-$$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$18 information from a scissored and stitched segment; the reported time cost is nearly independent of the target temperature, including extremely low temperatures (Wang et al., 2022). In open-system settings, T-TEDOPA achieves finite-temperature dynamics without explicit ancillas in the bosonic bath by absorbing temperature into a modified spectral density, while fermionic environments use a thermofield-like two-chain mapping (Lacroix et al., 2024).

A final source of ambiguity is terminological. In condensed-matter and statistical-physics usage, finite-temperature tensor network methods refer to approximations of thermal states, partition functions, or finite-$$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$19 open-system dynamics. In supervised learning, “finite temperature tensor network” denotes an MPS classifier modified by a Boltzmann-weighted temperature layer and a database-specific temperature parameter, not a solver for $$\langle O\rangle_\beta=\mathrm{Tr}(Oe^{-\beta H})/Z$$20 of a quantum lattice Hamiltonian (Lin et al., 2021). The shared language is the Boltzmann factor; the operational object is different.

Taken together, these developments indicate a field organized less by a single canonical formalism than by a recurring design problem: how to encode finite-temperature structure in a tensor network while keeping the relevant contractions, environments, and variational manifolds tractable. The direct-operator, purification, vectorization, renormalization, tailoring, and chain-mapping branches are best understood as complementary resolutions of that problem rather than mutually exclusive definitions of the method.