CFT spectrum for fermions

docs/src/cft.md

@@ -1,11 +1,12 @@
 # Conformal Field Theory Data
+
 TNRKit provides extensive tools for calculating conformal field theory data. Details about the implementation can be found in the TNRKit paper ([arxiv/2604.06922](https://arxiv.org/abs/2604.06922)).
 
-The core idea behind calculating the central charge, scaling dimensions, and conformal spins, is to calculate the spectrum of the fixed point tensor on a tube geometry. There are different ways to put fixed point tensors on a tube and the geometry of this tube is characterised by 3 parameters:
+The core idea behind calculating the central charge, scaling dimensions, and conformal spins is to calculate spectra of **transfer matrices** constructed from fixed point tensors on a tube geometry. There are different ways to put fixed point tensors on a tube, and the geometry of this tube is characterised by 3 parameters:
 $$[h, L, x]$$
-Where $h$ is the height of the tube, $L$ is the circumference, and $x$ is the horizontal shift. The higher the ratio $\frac{L}{h}$, the higher the resolution but also the more expensive the calculation.
+where $h$ is the height of the tube, $L$ is the circumference, and $x$ is the horizontal shift. The higher the ratio $\frac{L}{h}$, the higher the resolution but also the more expensive the calculation.
 
-To calculate cft data we provide the `CFTData` struct which can be used in the following ways:
+To calculate CFT data we provide the `CFTData` struct, which can be used in the following ways:
 
 ```julia
 CFTData(scheme; shape=[h, L, x])
@@ -15,18 +16,35 @@ CFTData(TA::TensorMap, TB::TensorMap; kwargs...) # 2x2 checkerboard unitcell
 
 The shapes we provide are: $[1, 1, 0]$, $[\sqrt{2}, 2\sqrt{2}, 0]$, $[1, 4, 1]$, $[1, 8, 1]$, $[\frac{4}{\sqrt{10}}, 2 \sqrt{10}, \frac{2}{\sqrt{10}}]$
 
-The last two of which require intermediate truncation steps, the parameters of which can be tuned by:
+The last two shapes require intermediate truncation steps, whose parameters can be tuned by:
+
 ```julia
 CFTData(scheme; shape=[1, 8, 1], trunc = trunc1, truncentanglement=trunc2)
 ```
 
-# CFTData struct
-The `CFTData` struct has two fields:
-- `central_charge`
-- `scaling_dimensions`
+## CFTData Struct
+
+The `CFTData` struct has three fields:
+
+- `central_charge`: a complex number which is the central charge of the CFT.
+- `modular_parameter`: a complex number which is the modular parameter of one tensor before building the tube transfer matrix.
+- `scaling_dimensions`: a `StructuredVector` storing the scaling dimensions and conformal spins of the primary fields and their descendants. It can be indexed like an `AbstractVector` (i.e. with scalars, slices, ...), or with sectors (e.g. `Z2Irrep(0)`), which will provide the scaling dimensions associated with that sector/charge.
+    To check which sectors you can index the `scaling_dimensions` with you can use `keys(scaling_dimensions)`.
+
+## How Fermionic CFT Data Is Extracted
 
-The `central_charge` can be either `missing` (when using the $[1, 1, 0]$ shape), or a number.
-The `scaling_dimensions` field is a `StructuredVector`.
+For a bosonic tensor network, `CFTData` builds a transfer matrix on the requested tube and diagonalizes it in each ordinary charge sector. For a fermionic tensor network, the same construction has one extra piece of global data: the **spin structure** around the tube. TNRKit extracts it by evaluating the same transfer matrix with two choices of boundary condition:
+
+- `:R`: periodic fermions around the tube, obtained with `pbc = true`. The macro `@tensor` automatically inserts a fermionic twist to take the *supertrace* across the tube.
+- `:NS`: antiperiodic fermions around the tube. Internally this is obtained by setting `pbc = false`, which explicitly inserts an additional fermionic twist to cancel the automatic supertrace twist, leaving the ordinary trace across the tube.
+
+The result is a `StructuredVector` whose keys are tuples `(spin_structure, charge)`. For fermionic systems without additional symmetries besides the fermion parity, there will be four sectors:
+
+```julia
+(:NS, FermionParity(0))
+(:NS, FermionParity(1))
+(:R,  FermionParity(0))
+(:R,  FermionParity(1))
+```
 
-The `scaling_dimensions` can be indexed like an `AbstractVector` (i.e. with scalars, slices, ...), or with sectors (e.g. `Z2Irrep(0)`), which will provide the scaling dimensions associated with that sector/charge.
-To check which sectors you can index the `scaling_dimensions` with you can use `keys(scaling_dimensions)`.
+The largest eigenvalue (corresponding to the identity field) should be in the NS even sector.

src/TNRKit.jl

@@ -15,6 +15,7 @@ using NonlinearSolve
 using Base.Threads
 using Combinatorics: permutations
 import TensorKitTensors.SpinOperators as SO
+import TensorKitTensors.FermionOperators as FO
 
 # stop criteria
 include("utility/stopping.jl")
@@ -125,7 +126,7 @@ export phi4_complex, phi4_complex_impϕ, phi4_complex_impϕdag, phi4_complex_imp
 
 include("models/quantum_1D.jl")
 export gate_to_tensor, vertical_stack_exp, vertical_stack_linear
-export quantum_ising_chain
+export quantum_ising_chain, kitaev_chain
 
 # utility functions
 include("utility/free_energy.jl")

src/models/quantum_1D.jl

@@ -74,6 +74,11 @@ end
 # =================================================
 
 """
+    quantum_ising_chain(dt::Float64; kwargs...)
+    quantum_ising_chain(elt::Type{<:Number}, dt::Float64; kwargs...)
+    quantum_ising_chain(symm::Type{<:Sector}, dt::Float64; kwargs...)
+    quantum_ising_chain(elt::Type{<:Number}, symm::Type{<:Sector}, dt::Float64; J::Float64=1.0, g::Float64=0.0)
+
 Partition function tensor for 1D transverse field Ising chain
 ```
     H(PBC) = -J ∑_i (σz_i σz_{i+1} + g σx_i)
@@ -95,3 +100,59 @@ quantum_ising_chain(elt::Type{<:Number}, dt::Float64; kwargs...) =
     quantum_ising_chain(elt, Trivial, dt; kwargs...)
 quantum_ising_chain(symm::Type{<:Sector}, dt::Float64; kwargs...) =
     quantum_ising_chain(ComplexF64, symm, dt; kwargs...)
+quantum_ising_chain(dt::Float64; kwargs...) =
+    quantum_ising_chain(ComplexF64, Trivial, dt; kwargs...)
+
+"""
+    kitaev_chain(dt::Float64; kwargs...)
+    kitaev_chain(elt::Type{<:Number}, dt::Float64; kwargs...)
+    kitaev_chain(symm::Type{<:Sector}, dt::Float64; kwargs...)
+    kitaev_chain(elt::Type{<:Number}, symm::Type{<:Sector}, dt::Float64; t::Float64=1.0, Δ::Float64=1.0, V::Float64=0.0, µ::Float64=0.0)
+
+Partition function tensor for 1D Kitaev chain model
+```
+    H = ∑_i [
+        (-t) (c†_i c_{i+1} + h.c.) + Δ (c_i c_{i+1} + h.c.)
+        + V (n_i - 1/2) (n_{i+1} - 1/2) - μ(n_i - 1/2)
+    ]
+```
+It is related to the spin-1/2 Heisenberg XYZ model
+```
+    H = ∑_i (J_x Sx_i Sx_{i+1} + J_y Sy_i Sy_{i+1}
+            + J_z Sz_i Sz_{i+1} - h Sz_j)
+```
+by Jordan-Wigner transformation
+```
+    t = -(Jx + Jy) / 4,  V = Jz,
+    Δ = -(Jx - Jy) / 4,  μ = h.
+```
+Special Cases
+- t = 1, Δ = 1, V = 0, μ = 2:   transverse field Ising model
+- t = 1, Δ = 0, V ≠ 0, μ = 0:   Heisenberg XXZ model
+"""
+function kitaev_chain(
+        elt::Type{<:Number}, symm::Type{<:Sector}, dt::Float64;
+        t::Float64 = 1.0, Δ::Float64 = 1.0, V::Float64 = 0.0, µ::Float64 = 0.0
+    )
+    fpfm = FO.f_plus_f_min(elt, symm)
+    hopping = (-t) * (fpfm + fpfm')
+    num = FO.f_num(elt, symm)
+    unit = TensorKit.id(codomain(num, 1))
+    num = num - 0.5 * unit
+    interac = V * (num ⊗ num)
+    chempot = -µ * (num ⊗ unit + unit ⊗ num) / 2
+    gate = hopping + interac + chempot
+    if Δ != 0
+        fmfm = FO.f_min_f_min(elt, symm)
+        pairing = Δ * (fmfm + fmfm')
+        gate = gate + pairing
+    end
+    gate = exp(-dt * gate)
+    return gate_to_tensor(gate)
+end
+kitaev_chain(elt::Type{<:Number}, dt::Float64; kwargs...) =
+    kitaev_chain(elt, Trivial, dt; kwargs...)
+kitaev_chain(symm::Type{<:Sector}, dt::Float64; kwargs...) =
+    kitaev_chain(ComplexF64, symm, dt; kwargs...)
+kitaev_chain(dt::Float64; kwargs...) =
+    kitaev_chain(ComplexF64, Trivial, dt; kwargs...)

src/schemes/looptnr.jl

@@ -9,12 +9,8 @@ Loop Optimization for Tensor Network Renormalization
     $(FUNCTIONNAME)(unitcell_2x2::Matrix{T})
 
 # Running the algorithm
-    run!(::LoopTNR, trunc::TruncationStrategy, criterion::stopcrit, parameters::LoopParameters, finalizer::Finalizer[,
-              entanglement_criterion::stopcrit, finalize_beginning=true, verbosity=1])
-
-    run!(::LoopTNR, trscheme::TruncationStrategy, criterion::stopcrit, parameters::LoopParameters; kwargs...)
-
-    run!(::LoopTNR, trscheme::TruncationStrategy, criterion::stopcrit[finalize_beginning=true, verbosity=1])
+    run!(::LoopTNR, trunc::TruncationStrategy, criterion::stopcrit, [parameters::LoopParameters], [finalizer::Finalizer];
+        [entanglement_criterion::stopcrit, finalize_beginning=true, verbosity=1])
 
 # LoopParameters
 See also: [`LoopParameters`](@ref)
@@ -490,7 +486,7 @@ function step!(
 end
 
 function run!(
-        scheme::LoopTNR, trscheme::TruncationStrategy,
+        scheme::LoopTNR, trunc::TruncationStrategy,
         criterion::stopcrit, loop_condition::LoopParameters,
         finalizer::Finalizer{E};
         entanglement_criterion = default_entanglement_criterion,
@@ -510,7 +506,7 @@ function run!(
 
         t = @elapsed while crit
             @infov 2 "Step $(steps + 1), data[end]: $(!isempty(data) ? data[end] : "empty")"
-            step!(scheme, trscheme, entanglement_criterion, loop_condition, verbosity)
+            step!(scheme, trunc, entanglement_criterion, loop_condition, verbosity)
             push!(data, finalizer.f!(scheme))
 
             steps += 1
@@ -522,21 +518,23 @@ function run!(
     return data
 end
 
-function run!(scheme, trscheme, criterion, loop_condition; kwargs...)
-    return run!(scheme, trscheme, criterion, loop_condition, default_Finalizer; kwargs...)
-end
-
-function run!(
-        scheme::LoopTNR, trscheme::TruncationStrategy, criterion::stopcrit;
-        finalize_beginning = true, verbosity = 1
-    )
-    loop_condition = LoopParameters()
-    return run!(
-        scheme, trscheme, criterion, loop_condition;
-        finalize_beginning = finalize_beginning,
-        verbosity = verbosity
-    )
-end
+# use default finalizer
+run!(
+    scheme::LoopTNR, trunc::TruncationStrategy,
+    criterion::stopcrit, loop_condition::LoopParameters; kwargs...
+) = run!(scheme, trunc, criterion, loop_condition, default_Finalizer; kwargs...)
+
+# use default loop parameters
+run!(
+    scheme::LoopTNR, trunc::TruncationStrategy,
+    criterion::stopcrit, finalizer::Finalizer; kwargs...
+) = run!(scheme, trunc, criterion, LoopParameters(), finalizer; kwargs...)
+
+# use both default loop paramater and default finalizer
+run!(
+    scheme::LoopTNR, trunc::TruncationStrategy,
+    criterion::stopcrit; kwargs...
+) = run!(scheme, trunc, criterion, LoopParameters(), default_Finalizer; kwargs...)
 
 function Base.show(io::IO, scheme::LoopTNR)
     println(io, "LoopTNR - Loop Tensor Network Renormalization")

src/utility/cft.jl

@@ -1,5 +1,5 @@
 """
-    struct CFTData{E, K, V, A <: AbstractVector{E}}
+    struct CFTData{E, K, A <: AbstractVector{E}}
 
 A struct to hold conformal data extracted from a TNR scheme.
 
@@ -11,16 +11,16 @@ A struct to hold conformal data extracted from a TNR scheme.
 # Fields
     - `central_charge::E`: The central charge of the CFT.
     - `modular_parameter::E`: The elementary modular parameter of a single tensor.
-    - `scaling_dimensions::StructuredVector{E, K, V, A}`: The scaling dimensions of the CFT, organized in a `StructuredVector` where the sectors correspond to different spin sectors (or other quantum numbers) and the data contains the scaling dimensions within those sectors
+    - `scaling_dimensions::StructuredVector{E, K, A}`: The scaling dimensions of the CFT, organized in a `StructuredVector` where the sectors correspond to different spin sectors (or other quantum numbers) and the data contains the scaling dimensions within those sectors
 
 """
-struct CFTData{E, K, V, A <: AbstractVector{E}}
+struct CFTData{E, K, A <: AbstractVector{E}}
     "Central charge of the CFT."
     central_charge::E
     "Elementary modular parameter for one tensor"
     modular_parameter::E
     "Scaling dimensions of the CFT."
-    scaling_dimensions::StructuredVector{E, K, V, A}
+    scaling_dimensions::StructuredVector{E, K, A}
 end
 
 function Base.show(io::IO, data::CFTData)
@@ -73,16 +73,7 @@ with `unitcell` copies of `T` concatenated horizontally.
 `τ0` is the modular parameter of a single `T`.
 """
 function _scaling_dimensions(T::TensorMap{E, S, 2, 2}, τ0::Number; unitcell = 1) where {E, S}
-    indices = [[i, -i, -(i + unitcell), i + 1] for i in 1:unitcell]
-    indices[end][4] = 1
-
-    T = ncon(fill(T, unitcell), indices)
-    # restore leg convention
-    outinds = Tuple(collect(1:unitcell))
-    ininds = Tuple(collect((unitcell + 1):(2unitcell)))
-    T = permute(T, (outinds, ininds))
-
-    sv = StructuredVector(eig_vals(T))
+    sv = _rowtm_eigvals(T, unitcell)
     sv = filter(x -> real(x) > 0 && abs(x) > 1.0e-12, sv)
     isempty(sv) && throw(ArgumentError("No valid eigenvalues found in transfer matrix spectrum."))
 
@@ -100,7 +91,8 @@ function area_term(
         TA::TensorMap{E, S, 2, 2}, TB::TensorMap{E, S, 2, 2}; is_real = true
     ) where {E, S}
     I = sectortype(TA)
-    λ = first(leading_eigenvalue(CFTTransferMatrix(TA, TB, [2, 2, 0]), one(I)))
+    pbc = (BraidingStyle(I) == Fermionic()) ? false : true
+    λ = first(leading_eigenvalue(CFTTransferMatrix(TA, TB, [2, 2, 0]), one(I); pbc, Nh = 1))
     return is_real ? real(λ) : λ
 end
 
@@ -120,18 +112,18 @@ function spec(
         τ0::Number; Nh = 25, trunc = notrunc(), truncentanglement = notrunc()
     ) where {E, S}
     I = sectortype(TA)
-    if BraidingStyle(I) != Bosonic()
-        throw(ArgumentError("Sectors with non-Bosonic charge $I has not been implemented"))
-    end
-
     tm = CFTTransferMatrix(TA, TB, shape; trunc, truncentanglement)
     τ = modular_parameter(tm, τ0)
 
     # eigenvalues of the transfer matrix from all charge sectors
     eigs = leading_eigenvalue(tm; Nh)
 
     # central charge
-    λ0 = eigs[one(I)][1]
+    λ0 = if BraidingStyle(I) == Fermionic()
+        eigs[(:NS, one(I))][1]
+    else
+        eigs[one(I)][1]
+    end
     area = shape[1] * shape[2]
     central_charge = 6 / pi / (imag(τ) - imag(τ0) * area / 4) * log(λ0)
 
@@ -158,8 +150,10 @@ end
 sigmoid(x) = 1 / (1 + exp(-x))
 logit(p) = log(p / (1 - p))
 function _find_λ0(TA, TB, shape)
-    charge = one(sectortype(TA))
-    λs = leading_eigenvalue(CFTTransferMatrix(TA, TB, shape), charge; Nh = 1)
+    I = sectortype(TA)
+    charge = one(I)
+    pbc = (BraidingStyle(I) == Fermionic()) ? false : true
+    λs = leading_eigenvalue(CFTTransferMatrix(TA, TB, shape), charge; pbc, Nh = 1)
     return real(first(λs))
 end
 

src/utility/structuredvector.jl

@@ -1,9 +1,9 @@
 """
-    StructuredVector{E, K, V, A} <: AbstractVector{E}
+    StructuredVector{E, K, A} <: AbstractVector{E}
 
 A vector whose elements are partitioned into named sectors.  Internally, data
-is stored as a flat `AbstractVector{E}` and a `Dict{K, V}` maps each sector key
-to the indices that belong to it.
+is stored as a flat `AbstractVector{E}` and a `Dict{K, Vector{Int}}` maps each
+sector key to the indices that belong to it.
 
 Supports the `AbstractVector` interface (integer indexing, `length`, `eachindex`,
 …), sector-based access via `v[sector]`, and the full `Dict` key interface
@@ -14,15 +14,15 @@ broadcasting with scalars all preserve the sector structure.
 
     StructuredVector(sv::SectorVector)
     StructuredVector(dict::Dict{K, <:AbstractVector{E}}) where {K, E}
-    StructuredVector(data::AbstractVector{E}, structure::Dict{K, V})
+    StructuredVector(data::AbstractVector{E}, structure::Dict{K, Vector{Int}})
 
 - From a TensorKit `SectorVector`.
 - From a dictionary mapping sectors to their data vectors.
 - Directly from a flat data array and a sector‑index mapping.
 """
-struct StructuredVector{E, K, V, A <: AbstractVector{E}} <: AbstractVector{E}
+struct StructuredVector{E, K, A <: AbstractVector{E}} <: AbstractVector{E}
     data::A
-    structure::Dict{K, V}
+    structure::Dict{K, Vector{Int}}
 end
 
 function StructuredVector(sv::TensorKit.SectorVector)
@@ -84,6 +84,35 @@ Base.:/(x::Number, v::StructuredVector) = StructuredVector(x ./ v.data, v.struct
 
 Base.keys(v::StructuredVector) = keys(v.structure)
 
+function Base.vcat(v1::StructuredVector{<:Any, K1}, v2::StructuredVector{<:Any, K2}) where {K1, K2}
+    new_data = vcat(v1.data, v2.data)
+    n1 = length(v1)
+    K = promote_type(K1, K2)
+    new_structure = Dict{K, Vector{Int}}()
+    for (k, inds) in v1.structure
+        new_structure[k] = copy(inds)
+    end
+    for (k, inds) in v2.structure
+        adjusted = inds .+ n1
+        if haskey(new_structure, k)
+            append!(new_structure[k], adjusted)
+        else
+            new_structure[k] = adjusted
+        end
+    end
+    return StructuredVector(new_data, new_structure)
+end
+
+# Julia Base provides a specialized `reduce(::typeof(vcat), A)` that bypasses
+# pairwise `vcat` calls and uses `similar` + bulk copy internally.  That path
+# does not know about sector structure and would produce a plain `Vector`.
+# We override it to fall back to pair-wise reduction which preserves the
+# StructuredVector container.
+function Base.reduce(::typeof(vcat), A::AbstractVector{<:StructuredVector})
+    isempty(A) && throw(ArgumentError("reducing vcat over an empty collection is not supported"))
+    return foldl(vcat, A)
+end
+
 # -- Broadcasting support -----------------------------------------------------
 # Custom broadcast style so that element-wise operations preserve the
 # StructuredVector container (and therefore the sector structure).
@@ -111,3 +140,15 @@ function Base.similar(bc::Broadcast.Broadcasted{StructuredVectorStyle}, ::Type{E
     end
     return StructuredVector(similar(sv.data, ElType), copy(sv.structure))
 end
+
+"""
+    mapkeys(f, v::StructuredVector)
+
+Return a new `StructuredVector` with each `structure` key transformed
+by a function `f`. The underlying `data` is shared with `v`.
+```
+"""
+function mapkeys(f, v::StructuredVector)
+    new_structure = Dict(f(k) => copy(inds) for (k, inds) in v.structure)
+    return StructuredVector(v.data, new_structure)
+end