[Feature discussion] Support for PBC FusionTree #594

Activity

1. lkdvos commented on Oct 26, 2026

   

   lkdvos

   on Oct 26, 2026

   Member

   More actions

   Hi,

   Thanks for checking out our packages.
   TensorKit indeed only handles canonical fusion trees, i.e. no PBC.
   The reason for this is that in general, the fusion diagrams you would draw with a closed loop forms an overcomplete set to expand a tensor in, while the canonical fusion trees form a complete orthonormal basis.
   This also means that you can always expand the "periodic trees" in terms of the canonical fusion trees, which is the recommended way to do this with TensorKit.

   For example, in order to do an ED calculation on 4 sites of a fibonacci-like chain, you can simply do the following:

   using TensorKit, LinearAlgebra
   
   N = 4                               # number of sites
   V = Vect[FibonacciAnyon](:τ => 1)   # one τ anyon per site
   𝟙 = id(ComplexF64, V)
   
   # Golden chain interaction: favour neighbouring τ's fusing to the vacuum
   h = zeros(ComplexF64, V ⊗ V ← V ⊗ V)
   block(h, FibonacciAnyon(:I)) .= -1
   
   # Periodic boundary conditions: the Hilbert space splits into sectors labelled by
   # the anyon flux threading the ring.
   
   ## Trivial flux sector
   # Ordinary open boundary tensors on V^N. For the boundary term h_{N,1} we move site N
   # to the front. In the trivial total charge sector, braiding over or under gives the same result.
   hᵢ(i) = reduce(⊗, [fill(𝟙, i - 1); h; fill(𝟙, N - 1 - i)])
   H = sum(hᵢ, 1:N-1)
   
   dom = ntuple(i -> N + i, N)
   T = braid(id(ComplexF64, V^N), ((N, 1:N-1...), dom), (1:N-1..., 100, dom...))
   H += T' * hᵢ(1) * T
   
   println("flux 1: ", eigvals(Hermitian(block(H, FibonacciAnyon(:I)))))
   
   # output
   flux 1: [-3.236067977499789, -0.7639320225002105]

   Things do get slightly more interesting if you also want to include "charged" states though.
   There, the typical trick is to insert an auxiliary index that carries this flux.
   For PBC, we then have to take care to correctly braid the periodic interaction around that leg though, for example the code below explicitly constructs the translation operator that braids around that flux leg:

   ## τ flux sector
   # Add an extra leg W that carries the flux. Physical states are those where V^N ⊗ W has total charge τ.
   W = Vect[FibonacciAnyon](:τ => 1)
   P = V^N ⊗ W
   
   hᵢ_τ(i) = hᵢ(i) ⊗ id(ComplexF64, W)
   H_τ = sum(hᵢ_τ, 1:N-1)
   
   # boundary term h_{N,1}: site N goes over the flux line, then back under it and the
   # other sites to the front, so it encircles the flux once.
   dom = ntuple(i -> N + 1 + i, N + 1)
   T = braid(id(ComplexF64, P), ((1:N-1..., N + 1, N), dom), (1:N-1..., 100, N, dom...))
   T = braid(T, ((N + 1, 1:N-1..., N), dom), (1:N..., -100, dom...))
   H_τ += T' * hᵢ_τ(1) * T
   
   println("flux τ: ", eigvals(Hermitian(block(H_τ, FibonacciAnyon(:τ)))))
   
   # output
   
   flux τ: [-2.999999999999999, -2.642684597778774, -1.0, -0.9999999999999993, -0.35731540222122404]

   ---

   If you are interested in going beyond ED, note that MPSKit.jl has full support for this:

   ## Same thing with DMRG, using MPSKit
   using MPSKit
   D = Vect[FibonacciAnyon](:I => 4, :τ => 4)   # maximal MPS bond space
   
   # trivial flux: the wrap-around term is just a long-range term between sites 1 and N
   lattice = fill(V, N)
   H_mpo = FiniteMPOHamiltonian(lattice, [(i, i + 1) => h for i in 1:N-1]..., (1, N) => h)
   ψ = FiniteMPS(lattice, D)
   ψ, = find_groundstate(ψ, H_mpo, DMRG())
   println("DMRG flux 1: ", real(expectation_value(ψ, H_mpo)))
   
   # τ flux: put the flux leg W as an extra site at the end, and fix the total charge
   # through the right virtual space of the MPS. The wrap-around term becomes a
   # three-site term on (1, N, N + 1): site N encircles the flux, then interacts with site 1.
   d = (4, 5, 6)
   M = braid(id(ComplexF64, V ⊗ V ⊗ W), ((1, 3, 2), d), (1, -100, 3, d...))
   M = braid(M, ((1, 3, 2), d), (1, 2, 100, d...))
   h_wrap = M' * (h ⊗ id(ComplexF64, W)) * M
   
   lattice_τ = [fill(V, N); W]
   H_mpo_τ = FiniteMPOHamiltonian(lattice_τ, [(i, i + 1) => h for i in 1:N-1]..., (1, N, N + 1) => h_wrap)
   ψ_τ = FiniteMPS(lattice_τ, D; right = W)
   ψ_τ, = find_groundstate(ψ_τ, H_mpo_τ, DMRG())
   println("DMRG flux τ: ", real(expectation_value(ψ_τ, H_mpo_τ)))
   
   # output
   DMRG flux 1: -3.2361
   DMRG flux τ: -3.0000