From 2b76bd33b4665e54e31777b4b3d24c76a5833ce4 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Mon, 8 Jun 2026 17:26:29 +0800
Subject: [PATCH 01/14] Add (elementary) modular parameter to CFTData

---
 src/utility/cft.jl | 118 ++++++++++++++++++++++++++++++++++++---------
 1 file changed, 96 insertions(+), 22 deletions(-)

diff --git a/src/utility/cft.jl b/src/utility/cft.jl
index ebcf768f..db668710 100644
--- a/src/utility/cft.jl
+++ b/src/utility/cft.jl
@@ -55,9 +55,10 @@ A struct to hold conformal data extracted from a TNR scheme.
 
 """
 struct CFTData{E, I}
-    "Central charge of the CFT. Will be `nothing` if not calculated."
+    "Central charge of the CFT. Will be `missing` if not calculated."
     central_charge::Union{E, Missing}
-
+    "Elementary modular parameter for one tensor"
+    modular_param::E
     "Scaling dimensions of the CFT."
     scaling_dimensions::TensorKit.SectorVector{E, I}
 end
@@ -77,7 +78,7 @@ function CFTData(T::TensorMap{E, S, 2, 2}; shape = [sqrt(2), 2 * sqrt(2), 0], kw
     end
 end
 CFTData(scheme::TNRScheme; kwargs...) = CFTData(scheme.T; kwargs...) # simple 1-site unitcell schemes
-CFTData(scheme::LoopTNR; kwargs...) = CFTData(scheme.TA, scheme.TB; kwargs...) # simple 1-site unitcell schemes
+CFTData(scheme::LoopTNR; kwargs...) = CFTData(scheme.TA, scheme.TB; kwargs...) # 2-site unitcell schemes
 function CFTData(scheme::BTRG; kwargs...) # merge bond tensors into central tensor
     @tensor T_unit[-1 -2; -3 -4] := scheme.T[1 2; -3 -4] * scheme.S1[-2; 2] *
         scheme.S2[-1; 1]
@@ -86,18 +87,14 @@ end
 
 # Main implementation, two-site unitcell
 function CFTData(TA::TensorMap{E, S, 2, 2}, TB::TensorMap{E, S, 2, 2}; shape = [sqrt(2), 2 * sqrt(2), 0], trunc = truncrank(16), truncentanglement = trunctol(; rtol = 1.0e-14)) where {E, S}
+    norm_const = area_term(TA, TB)^(1 / 4) # canonical normalisation constant
+    TA′, TB′ = TA / norm_const, TB / norm_const
     if shape == [1, 1, 0]
         throw(ArgumentError("The shape [1, 1, 0] is not compatible with a two-site unit cell."))
     elseif (shape ≈ [sqrt(2), 2 * sqrt(2), 0]) || (shape == [1, 4, 1]) # these shapes need no truncation
-        norm_const = area_term(TA, TB)^(1 / 4) # canonical normalisation constant
-        return spec(TA / norm_const, TB / norm_const, shape)
+        return spec(TA′, TB′, shape)
     elseif (shape == [1, 8, 1]) || (shape ≈ [4 / sqrt(10), 2 * sqrt(10), 2 / sqrt(10)])
-
-        norm_const = area_term(TA, TB)^(1 / 4) # canonical normalisation constant
-
-        dl, ur, ul, dr = MPO_opt(
-            TA / norm_const, TB / norm_const, trunc, truncentanglement
-        )
+        dl, ur, ul, dr = MPO_opt(TA′, TB′, trunc, truncentanglement)
         T = reduced_MPO(dl, ur, ul, dr, trunc)
         return spec(T, T, shape)
     else
@@ -129,20 +126,16 @@ end
 
 """
 The "canonical" normalization constant for loop-TNR tensors,
-which is the eigenvalue with largest real part of the 2 x 2 transfer matrix.
+which is the eigenvalue with largest norm of the 2 x 2 transfer matrix.
 """
 function area_term(A, B; is_real = true)
-    a_in = domain(A)[1]
-    b_in = domain(B)[1]
-    x0 = ones(a_in ⊗ b_in)
-
     function f0(x)
         @plansor fx[-1 -2] := A[c -1; 1 m] * x[1 2] * B[m -2; 2 c]
         @plansor ffx[-1 -2] := B[c -1; 1 m] * fx[1 2] * A[m -2; 2 c]
         return ffx
     end
-
-    spec0, _, info = eigsolve(f0, x0, 1, :LR; verbosity = 0)
+    x0 = ones(domain(A, 1) ⊗ domain(B, 1))
+    spec0, _, info = eigsolve(f0, x0, 1, :LM; verbosity = 0)
     if info.converged == 0
         @warn "The area term eigensolver did not converge."
     end
@@ -155,9 +148,10 @@ end
 
 # The case with spin is based on https://arxiv.org/pdf/1512.03846 and some private communications with Yingjie Wei and Atsushi Ueda
 function spec(TA::TensorMap, TB::TensorMap, shape::Array; Nh = 25)
-    area = shape[1] * shape[2]
-    Imτ = shape[1] / shape[2]
-    relative_shift = shape[3] / shape[1]
+    τ0 = _modular_parameter(TA, TB)
+    area = imag(τ0) * shape[1] * shape[2]
+    Imτ = imag(τ0) * shape[1] / shape[2]
+    relative_shift = real(τ0) / imag(τ0) + shape[3] / (shape[1] * imag(τ0))
 
     I = sectortype(TA)
     𝔽 = field(TA)
@@ -219,7 +213,7 @@ function spec(TA::TensorMap, TB::TensorMap, shape::Array; Nh = 25)
     sv = sort(sv; by = real)
     sv = filter(x -> real(x) ≤ 1.0e16, sv)
 
-    return CFTData(central_charge, sv)
+    return CFTData(central_charge, τ0, sv)
 end
 
 function MPO_opt(
@@ -281,3 +275,83 @@ function reduced_MPO(
     @plansor translate[-1 -2; -3 -4] := U[-2; 1 -4] * D[-1 1; -3]
     return translate
 end
+
+# Elementary modular parameter
+# ============================
+"""
+    _modular_parameter(TA::TensorMap{E, S, 2, 2}, TB::TensorMap{E, S, 2, 2})
+
+Extract the elementary modular parameter of one tensor from 2x2 transfer matrices.
+"""
+function _modular_parameter(TA::TensorMap{E, S, 2, 2}, TB::TensorMap{E, S, 2, 2}) where {E, S}
+    # vertical (north to south)
+    evs_v = _eigsolve_2x2_NtoS(TA, TB)
+    # horizontal (east to west)
+    evs_h = _eigsolve_2x2_EtoW(TA, TB)
+    # diagonal (northeast to southwest)
+    evs_a = _eigsolve_2x2_NEtoSW(TA, TB)
+    # diagonal (northwest to southeast)
+    evs_b = _eigsolve_2x2_NWtoSE(TA, TB)
+    # norm²
+    v = sqrt(
+        log(abs(evs_v[2] / evs_v[1])) /
+            log(abs(evs_h[2] / evs_h[1]))
+    )
+    # phase
+    r = log(abs(evs_a[2] / evs_a[1])) /
+        log(abs(evs_b[2] / evs_b[1]))
+    θ = acos((v^2 + 1) / (2 * v) * (r - 1) / (r + 1))
+    return v * cis(θ)
+end
+
+function _eigsolve_2x2_NtoS(TA, TB)
+    function f0(x)
+        @plansor fx[-1 -2] := TA[c -1; 1 m] * x[1 2] * TB[m -2; 2 c]
+        @plansor fx[-1 -2] := TB[c -1; 1 m] * fx[1 2] * TA[m -2; 2 c]
+        return fx
+    end
+    x0 = ones(domain(TA, 1) ⊗ domain(TB, 1))
+    spec0, _, info = eigsolve(f0, x0, 2, :LM; verbosity = 0)
+    if info.converged == 0
+        @warn "The area term eigensolver did not converge."
+    end
+    return spec0
+end
+
+function _eigsolve_2x2_EtoW(TA, TB)
+    TA′ = permute(TA, ((3, 1), (4, 2)))
+    TB′ = permute(TB, ((3, 1), (4, 2)))
+    return _eigsolve_2x2_NtoS(TA′, TB′)
+end
+
+function _eigsolve_2x2_NEtoSW(TA, TB)
+    #= 
+        1   2
+        ┌---┬---┐
+        |   |   |
+    3'--A---B---┤ 3
+        |   |   |
+    4'--B---A---┘ 4
+        |   |
+        1'  2'
+    =#
+    function f0(x)
+        @plansor fx[-1 -2 -3 -4] := TB[-2 -3; a b] * x[-1 a b -4]
+        @plansor fx[-1 -2 -3 -4] := TA[-1 -2; a b] * fx[a b -3 -4]
+        @plansor fx[-1 -2 -3 -4] := TA[-3 -4; a b] * fx[-1 -2 a b]
+        @plansor fx[-1 -2 -3 -4] := TB[-4 -1; a b] * fx[-3 a b -2]
+        return fx
+    end
+    x0 = ones(domain(TA, 1) ⊗ domain(TB, 1) ⊗ domain(TB, 2) ⊗ domain(TA, 2))
+    spec0, _, info = eigsolve(f0, x0, 2, :LM; verbosity = 0)
+    if info.converged == 0
+        @warn "The area term eigensolver did not converge."
+    end
+    return spec0
+end
+
+function _eigsolve_2x2_NWtoSE(TA, TB)
+    TA′ = permute(TA, ((2, 4), (1, 3)))
+    TB′ = permute(TB, ((2, 4), (1, 3)))
+    return _eigsolve_2x2_NEtoSW(TA′, TB′)
+end

From cef20e884566e0922184fa485b74b94c2154f569 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Tue, 9 Jun 2026 09:37:59 +0800
Subject: [PATCH 02/14] Comment on how lattice is rotated in TRG and loop-TNR

---
 src/schemes/looptnr.jl | 3 ++-
 src/schemes/trg.jl     | 3 ++-
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/schemes/looptnr.jl b/src/schemes/looptnr.jl
index e1d1140b..c077c487 100644
--- a/src/schemes/looptnr.jl
+++ b/src/schemes/looptnr.jl
@@ -433,7 +433,8 @@ function loop_opt(
 end
 
 """
-Coarse-grain `ΨB` to renormalized `TA`, `TB` tensors
+Coarse-grain `ΨB` to renormalized `TA`, `TB` tensors.
+The lattice is rotated by 135 degrees in counter clockwise direction.
 """
 function ΨB_to_TATB(psiB::Vector{T}) where {T <: AbstractTensorMap{<:Any, <:Any, 1, 2}}
     #= 
diff --git a/src/schemes/trg.jl b/src/schemes/trg.jl
index 021379e4..24a585e1 100644
--- a/src/schemes/trg.jl
+++ b/src/schemes/trg.jl
@@ -9,7 +9,8 @@ Tensor Renormalization Group
 # Running the algorithm
     run!(::TRG, trunc::TruncationStrategy, stop::Stopcrit[, finalizer=default_Finalizer, finalize_beginning=true, verbosity=1])
 
-Each step rescales the lattice by a (linear) factor of √2
+Each step rescales the lattice by a (linear) factor of √2,
+and rotate the lattice by 45 degrees in counter clockwise direction.
 
 !!! info "verbosity levels"
     - 0: No output

From 917f1aa75eec459df8954b532f42f9d34e29d25c Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Tue, 9 Jun 2026 12:06:10 +0800
Subject: [PATCH 03/14] Refactors

---
 src/utility/cft.jl | 35 ++++++++++++++++++++---------------
 1 file changed, 20 insertions(+), 15 deletions(-)

diff --git a/src/utility/cft.jl b/src/utility/cft.jl
index db668710..a0b091bd 100644
--- a/src/utility/cft.jl
+++ b/src/utility/cft.jl
@@ -72,7 +72,8 @@ end
 
 function CFTData(T::TensorMap{E, S, 2, 2}; shape = [sqrt(2), 2 * sqrt(2), 0], kwargs...) where {E, S}
     if shape == [1, 1, 0] # trivial implementation
-        return CFTData(missing, _scaling_dimensions(T))
+        τ0 = elementary_modular_parameter(T, T)
+        return CFTData(missing, τ0, _scaling_dimensions(T))
     else
         CFTData(T, T; shape, kwargs...)
     end
@@ -102,8 +103,7 @@ function CFTData(TA::TensorMap{E, S, 2, 2}, TB::TensorMap{E, S, 2, 2}; shape = [
     end
 end
 
-# Trivial diagonalisation of the transfer matrix. Currently the v and unitcell are not acessible from the outside.
-# The user should really be using the other shapes anyways.
+# TODO: replace v with solved elementary modular parameter
 function _scaling_dimensions(T::TensorMap{E, S, 2, 2}; v = 1, unitcell = 1) where {E, S}
     # stack unitcell copies of T and trace
     indices = [[i, -i, -(i + unitcell), i + 1] for i in 1:unitcell]
@@ -128,13 +128,13 @@ end
 The "canonical" normalization constant for loop-TNR tensors,
 which is the eigenvalue with largest norm of the 2 x 2 transfer matrix.
 """
-function area_term(A, B; is_real = true)
+function area_term(TA, TB; is_real = true)
     function f0(x)
-        @plansor fx[-1 -2] := A[c -1; 1 m] * x[1 2] * B[m -2; 2 c]
-        @plansor ffx[-1 -2] := B[c -1; 1 m] * fx[1 2] * A[m -2; 2 c]
+        @plansor fx[-1 -2] := TA[c -1; 1 m] * x[1 2] * TB[m -2; 2 c]
+        @plansor ffx[-1 -2] := TB[c -1; 1 m] * fx[1 2] * TA[m -2; 2 c]
         return ffx
     end
-    x0 = ones(domain(A, 1) ⊗ domain(B, 1))
+    x0 = ones(domain(TA, 1) ⊗ domain(TB, 1))
     spec0, _, info = eigsolve(f0, x0, 1, :LM; verbosity = 0)
     if info.converged == 0
         @warn "The area term eigensolver did not converge."
@@ -148,7 +148,7 @@ end
 
 # The case with spin is based on https://arxiv.org/pdf/1512.03846 and some private communications with Yingjie Wei and Atsushi Ueda
 function spec(TA::TensorMap, TB::TensorMap, shape::Array; Nh = 25)
-    τ0 = _modular_parameter(TA, TB)
+    τ0 = elementary_modular_parameter(TA, TB)
     area = imag(τ0) * shape[1] * shape[2]
     Imτ = imag(τ0) * shape[1] / shape[2]
     relative_shift = real(τ0) / imag(τ0) + shape[3] / (shape[1] * imag(τ0))
@@ -159,6 +159,7 @@ function spec(TA::TensorMap, TB::TensorMap, shape::Array; Nh = 25)
         throw(ArgumentError("Sectors with non-Bosonic charge $I has not been implemented"))
     end
 
+    # eigenvalues of the transfer matrix
     xspace, f = if shape ≈ [1, 4, 1]
         domain(TA)[1] ⊗ domain(TB)[1] ⊗ domain(TA)[1] ⊗ domain(TB)[1],
             MPO_action_1x4_twist
@@ -169,7 +170,6 @@ function spec(TA::TensorMap, TB::TensorMap, shape::Array; Nh = 25)
             shape ≈ [4 / sqrt(10), 2 * sqrt(10), 2 / sqrt(10)]
         domain(TB) ⊗ domain(TB), MPO_action_2gates
     end
-
     spec_sector = Dict(
         map(sectors(fuse(xspace))) do charge
             V = (I == Trivial) ? 𝔽^1 : Vect[I](charge => 1)
@@ -190,19 +190,21 @@ function spec(TA::TensorMap, TB::TensorMap, shape::Array; Nh = 25)
         end
     )
 
+    # central charge
     norm_const_0 = spec_sector[one(I)][1]
     central_charge = 6 / pi / (Imτ - area / 4) * log(norm_const_0)
 
-    # Construct a SectorVector from the data of the different sectors
+    # A SectorVector for scaling dimension and conformal spin
     data = ComplexF64[]
     structure = TensorKit.SectorDict{sectortype(xspace), UnitRange{Int}}()
     last_index = 1
     for charge in sectors(fuse(xspace))
         DeltaS = -1 / (2 * pi * Imτ) * log.(spec_sector[charge] / norm_const_0)
-        if !(relative_shift ≈ 0)
+        if !isapprox(relative_shift, 0; atol = 1.0e-6)
             push!(data, (real.(DeltaS) + imag.(DeltaS) / relative_shift * im)...)
             structure[charge] = last_index:(last_index + length(DeltaS) - 1)
         else
+            # not enough precision to resolve conformal spin
             push!(data, real.(DeltaS)...)
             structure[charge] = last_index:(last_index + length(DeltaS) - 1)
         end
@@ -278,12 +280,15 @@ end
 
 # Elementary modular parameter
 # ============================
+# TODO: one-tensor version
 """
-    _modular_parameter(TA::TensorMap{E, S, 2, 2}, TB::TensorMap{E, S, 2, 2})
+    elementary_modular_parameter(TA::TensorMap{E, S, 2, 2}, TB::TensorMap{E, S, 2, 2}) where {E, S}
 
 Extract the elementary modular parameter of one tensor from 2x2 transfer matrices.
 """
-function _modular_parameter(TA::TensorMap{E, S, 2, 2}, TB::TensorMap{E, S, 2, 2}) where {E, S}
+function elementary_modular_parameter(
+        TA::TensorMap{E, S, 2, 2}, TB::TensorMap{E, S, 2, 2}
+    ) where {E, S}
     # vertical (north to south)
     evs_v = _eigsolve_2x2_NtoS(TA, TB)
     # horizontal (east to west)
@@ -321,7 +326,7 @@ end
 function _eigsolve_2x2_EtoW(TA, TB)
     TA′ = permute(TA, ((3, 1), (4, 2)))
     TB′ = permute(TB, ((3, 1), (4, 2)))
-    return _eigsolve_2x2_NtoS(TA′, TB′)
+    return _eigsolve_2x2_NtoS(TB′, TA′)
 end
 
 function _eigsolve_2x2_NEtoSW(TA, TB)
@@ -353,5 +358,5 @@ end
 function _eigsolve_2x2_NWtoSE(TA, TB)
     TA′ = permute(TA, ((2, 4), (1, 3)))
     TB′ = permute(TB, ((2, 4), (1, 3)))
-    return _eigsolve_2x2_NEtoSW(TA′, TB′)
+    return _eigsolve_2x2_NEtoSW(TB′, TA′)
 end

From bd69abb2908e27486559362f4533ded9a6324941 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Tue, 9 Jun 2026 14:18:49 +0800
Subject: [PATCH 04/14] =?UTF-8?q?Note=20on=20how=20=CF=840=20changed=20in?=
 =?UTF-8?q?=20TRG=20and=20loop-TNR?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/schemes/looptnr.jl | 1 +
 src/schemes/trg.jl     | 1 +
 2 files changed, 2 insertions(+)

diff --git a/src/schemes/looptnr.jl b/src/schemes/looptnr.jl
index c077c487..8c5db9f1 100644
--- a/src/schemes/looptnr.jl
+++ b/src/schemes/looptnr.jl
@@ -435,6 +435,7 @@ end
 """
 Coarse-grain `ΨB` to renormalized `TA`, `TB` tensors.
 The lattice is rotated by 135 degrees in counter clockwise direction.
+The elementary modular parameter `τ₀ ↦ (1 + τ₀) / (1 - τ₀)`.
 """
 function ΨB_to_TATB(psiB::Vector{T}) where {T <: AbstractTensorMap{<:Any, <:Any, 1, 2}}
     #= 
diff --git a/src/schemes/trg.jl b/src/schemes/trg.jl
index 24a585e1..af1c9f43 100644
--- a/src/schemes/trg.jl
+++ b/src/schemes/trg.jl
@@ -11,6 +11,7 @@ Tensor Renormalization Group
 
 Each step rescales the lattice by a (linear) factor of √2,
 and rotate the lattice by 45 degrees in counter clockwise direction.
+The elementary modular parameter `τ₀ ↦ (τ₀ - 1) / (τ₀ + 1)`.
 
 !!! info "verbosity levels"
     - 0: No output

From 38e542ac8b2f353940a4200c478f35fc7559af2e Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Tue, 9 Jun 2026 15:09:31 +0800
Subject: [PATCH 05/14] Use full name "modular_parameter" in CFTData

---
 src/utility/cft.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/utility/cft.jl b/src/utility/cft.jl
index a0b091bd..aaa82578 100644
--- a/src/utility/cft.jl
+++ b/src/utility/cft.jl
@@ -51,6 +51,7 @@ A struct to hold conformal data extracted from a TNR scheme.
 
 # Fields
     - `central_charge::Union{E, Missing}`: The central charge of the CFT. Will be `nothing` if not calculated.
+    - `modular_parameter::E`: The elementary modular parameter of a square spacetime patch of the CFT.
     - `scaling_dimensions::TensorKit.SectorVector{E, I}`: The scaling dimensions of the CFT, organized in a `TensorKit.SectorVector` where the sectors correspond to different spin sectors (or other quantum numbers) and the data contains the scaling dimensions within those sectors
 
 """
@@ -58,7 +59,7 @@ struct CFTData{E, I}
     "Central charge of the CFT. Will be `missing` if not calculated."
     central_charge::Union{E, Missing}
     "Elementary modular parameter for one tensor"
-    modular_param::E
+    modular_parameter::E
     "Scaling dimensions of the CFT."
     scaling_dimensions::TensorKit.SectorVector{E, I}
 end

From cb51116682c9c1fb3ce3eb295532b3c7e8b548fb Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Tue, 9 Jun 2026 17:21:33 +0800
Subject: [PATCH 06/14] Fix for [sqrt(2), 2 * sqrt(2), 0]

---
 src/utility/cft.jl | 30 +++++++++++++++++-------------
 1 file changed, 17 insertions(+), 13 deletions(-)

diff --git a/src/utility/cft.jl b/src/utility/cft.jl
index aaa82578..3109e963 100644
--- a/src/utility/cft.jl
+++ b/src/utility/cft.jl
@@ -149,11 +149,6 @@ end
 
 # The case with spin is based on https://arxiv.org/pdf/1512.03846 and some private communications with Yingjie Wei and Atsushi Ueda
 function spec(TA::TensorMap, TB::TensorMap, shape::Array; Nh = 25)
-    τ0 = elementary_modular_parameter(TA, TB)
-    area = imag(τ0) * shape[1] * shape[2]
-    Imτ = imag(τ0) * shape[1] / shape[2]
-    relative_shift = real(τ0) / imag(τ0) + shape[3] / (shape[1] * imag(τ0))
-
     I = sectortype(TA)
     𝔽 = field(TA)
     if BraidingStyle(I) != Bosonic()
@@ -161,15 +156,20 @@ function spec(TA::TensorMap, TB::TensorMap, shape::Array; Nh = 25)
     end
 
     # eigenvalues of the transfer matrix
-    xspace, f = if shape ≈ [1, 4, 1]
+    τ0 = elementary_modular_parameter(TA, TB)
+    xspace, f, τ = if shape ≈ [1, 4, 1]
         domain(TA)[1] ⊗ domain(TB)[1] ⊗ domain(TA)[1] ⊗ domain(TB)[1],
-            MPO_action_1x4_twist
+            MPO_action_1x4_twist, (1 + τ0) / 4
     elseif shape ≈ [1, 8, 1]
         domain(TA)[1] ⊗ domain(TB)[1] ⊗ domain(TA)[1] ⊗ domain(TB)[1],
-            MPO_action_1x4
-    elseif shape ≈ [sqrt(2), 2 * sqrt(2), 0] ||
-            shape ≈ [4 / sqrt(10), 2 * sqrt(10), 2 / sqrt(10)]
-        domain(TB) ⊗ domain(TB), MPO_action_2gates
+            MPO_action_1x4, (1 + τ0) / 8
+    elseif shape ≈ [sqrt(2), 2 * sqrt(2), 0]
+        domain(TB) ⊗ domain(TB), MPO_action_2gates, (1 + τ0) / 2 / (1 - τ0)
+    elseif shape ≈ [4 / sqrt(10), 2 * sqrt(10), 2 / sqrt(10)]
+        error("Not implemented")
+        # domain(TB) ⊗ domain(TB), MPO_action_2gates
+    else
+        error("Unsupported transfer matrix shape.")
     end
     spec_sector = Dict(
         map(sectors(fuse(xspace))) do charge
@@ -193,15 +193,19 @@ function spec(TA::TensorMap, TB::TensorMap, shape::Array; Nh = 25)
 
     # central charge
     norm_const_0 = spec_sector[one(I)][1]
-    central_charge = 6 / pi / (Imτ - area / 4) * log(norm_const_0)
+    area = shape[1] * shape[2]
+    central_charge = 6 / pi / (imag(τ) - imag(τ0) * area / 4) * log(norm_const_0)
 
     # A SectorVector for scaling dimension and conformal spin
     data = ComplexF64[]
     structure = TensorKit.SectorDict{sectortype(xspace), UnitRange{Int}}()
     last_index = 1
+    relative_shift = real(τ) / imag(τ)
     for charge in sectors(fuse(xspace))
-        DeltaS = -1 / (2 * pi * Imτ) * log.(spec_sector[charge] / norm_const_0)
+        # DeltaS = Δ - i s Re(τ) / Im(τ)
+        DeltaS = -1 / (2 * pi * imag(τ)) * log.(spec_sector[charge] / norm_const_0)
         if !isapprox(relative_shift, 0; atol = 1.0e-6)
+            # save `Δ - i s` in `data`
             push!(data, (real.(DeltaS) + imag.(DeltaS) / relative_shift * im)...)
             structure[charge] = last_index:(last_index + length(DeltaS) - 1)
         else

From 48133c9b3c78d68db3951c4ff00642f403469ee1 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Tue, 9 Jun 2026 20:35:56 +0800
Subject: [PATCH 07/14] Fix for [1, 8, 1] and [4 / sqrt(10), 2 * sqrt(10), 2 /
 sqrt(10)]

---
 src/utility/cft.jl | 74 +++++++++++++++++++++++++++++++++-------------
 1 file changed, 54 insertions(+), 20 deletions(-)

diff --git a/src/utility/cft.jl b/src/utility/cft.jl
index 3109e963..8efa07d3 100644
--- a/src/utility/cft.jl
+++ b/src/utility/cft.jl
@@ -71,6 +71,15 @@ function Base.show(io::IO, data::CFTData)
     return nothing
 end
 
+CFTData(scheme::TNRScheme; kwargs...) = CFTData(scheme.T; kwargs...) # simple 1-site unitcell schemes
+CFTData(scheme::LoopTNR; kwargs...) = CFTData(scheme.TA, scheme.TB; kwargs...) # 2-site unitcell schemes
+function CFTData(scheme::BTRG; kwargs...) # merge bond tensors into central tensor
+    @tensor T_unit[-1 -2; -3 -4] := scheme.T[1 2; -3 -4] * scheme.S1[-2; 2] *
+        scheme.S2[-1; 1]
+    return CFTData(T_unit; kwargs...)
+end
+
+# one-site unitcell
 function CFTData(T::TensorMap{E, S, 2, 2}; shape = [sqrt(2), 2 * sqrt(2), 0], kwargs...) where {E, S}
     if shape == [1, 1, 0] # trivial implementation
         τ0 = elementary_modular_parameter(T, T)
@@ -79,26 +88,23 @@ function CFTData(T::TensorMap{E, S, 2, 2}; shape = [sqrt(2), 2 * sqrt(2), 0], kw
         CFTData(T, T; shape, kwargs...)
     end
 end
-CFTData(scheme::TNRScheme; kwargs...) = CFTData(scheme.T; kwargs...) # simple 1-site unitcell schemes
-CFTData(scheme::LoopTNR; kwargs...) = CFTData(scheme.TA, scheme.TB; kwargs...) # 2-site unitcell schemes
-function CFTData(scheme::BTRG; kwargs...) # merge bond tensors into central tensor
-    @tensor T_unit[-1 -2; -3 -4] := scheme.T[1 2; -3 -4] * scheme.S1[-2; 2] *
-        scheme.S2[-1; 1]
-    return CFTData(T_unit; kwargs...)
-end
 
 # Main implementation, two-site unitcell
-function CFTData(TA::TensorMap{E, S, 2, 2}, TB::TensorMap{E, S, 2, 2}; shape = [sqrt(2), 2 * sqrt(2), 0], trunc = truncrank(16), truncentanglement = trunctol(; rtol = 1.0e-14)) where {E, S}
+function CFTData(
+        TA::TensorMap{E, S, 2, 2}, TB::TensorMap{E, S, 2, 2}; shape = [sqrt(2), 2 * sqrt(2), 0],
+        trunc = truncrank(16), truncentanglement = trunctol(; rtol = 1.0e-14)
+    ) where {E, S}
     norm_const = area_term(TA, TB)^(1 / 4) # canonical normalisation constant
     TA′, TB′ = TA / norm_const, TB / norm_const
+    τ0 = elementary_modular_parameter(TA′, TB′)
     if shape == [1, 1, 0]
         throw(ArgumentError("The shape [1, 1, 0] is not compatible with a two-site unit cell."))
     elseif (shape ≈ [sqrt(2), 2 * sqrt(2), 0]) || (shape == [1, 4, 1]) # these shapes need no truncation
-        return spec(TA′, TB′, shape)
+        return spec(TA′, TB′, shape, τ0)
     elseif (shape == [1, 8, 1]) || (shape ≈ [4 / sqrt(10), 2 * sqrt(10), 2 / sqrt(10)])
         dl, ur, ul, dr = MPO_opt(TA′, TB′, trunc, truncentanglement)
         T = reduced_MPO(dl, ur, ul, dr, trunc)
-        return spec(T, T, shape)
+        return spec(T, T, shape, τ0)
     else
         throw(ArgumentError("Shape $shape is not implemented."))
     end
@@ -148,7 +154,19 @@ function area_term(TA, TB; is_real = true)
 end
 
 # The case with spin is based on https://arxiv.org/pdf/1512.03846 and some private communications with Yingjie Wei and Atsushi Ueda
-function spec(TA::TensorMap, TB::TensorMap, shape::Array; Nh = 25)
+"""
+    spec(TA::TensorMap, TB::TensorMap, shape::Array; Nh = 25)
+
+Internal function to construct transfer matrices and extract conformal data.
+
+# Parameters
+
+TA, TB: Tensors used to construct the transfer matrix. 
+    They can be different from the original tensor in the network.
+shape: A triplet `[h, L, x]` describing the geometry of the transfer matrix.
+τ0: The modular parameter of a square patch of the original network.
+"""
+function spec(TA::TensorMap, TB::TensorMap, shape::Array, τ0::Number; Nh = 25)
     I = sectortype(TA)
     𝔽 = field(TA)
     if BraidingStyle(I) != Bosonic()
@@ -156,18 +174,19 @@ function spec(TA::TensorMap, TB::TensorMap, shape::Array; Nh = 25)
     end
 
     # eigenvalues of the transfer matrix
-    τ0 = elementary_modular_parameter(TA, TB)
     xspace, f, τ = if shape ≈ [1, 4, 1]
         domain(TA)[1] ⊗ domain(TB)[1] ⊗ domain(TA)[1] ⊗ domain(TB)[1],
             MPO_action_1x4_twist, (1 + τ0) / 4
-    elseif shape ≈ [1, 8, 1]
-        domain(TA)[1] ⊗ domain(TB)[1] ⊗ domain(TA)[1] ⊗ domain(TB)[1],
-            MPO_action_1x4, (1 + τ0) / 8
     elseif shape ≈ [sqrt(2), 2 * sqrt(2), 0]
         domain(TB) ⊗ domain(TB), MPO_action_2gates, (1 + τ0) / 2 / (1 - τ0)
+        # in the following cases, (TA, TB) are no longer the original tensor in the network
+    elseif shape ≈ [1, 8, 1]
+        τ0′ = (τ0 - 1) / 2
+        domain(TA)[1] ⊗ domain(TB)[1] ⊗ domain(TA)[1] ⊗ domain(TB)[1],
+            MPO_action_1x4, τ0′ / 4
     elseif shape ≈ [4 / sqrt(10), 2 * sqrt(10), 2 / sqrt(10)]
-        error("Not implemented")
-        # domain(TB) ⊗ domain(TB), MPO_action_2gates
+        τ0′ = (τ0 - 1) / 2
+        domain(TB) ⊗ domain(TB), MPO_action_2gates, (1 + τ0′) / 2 / (1 - τ0′)
     else
         error("Unsupported transfer matrix shape.")
     end
@@ -252,6 +271,10 @@ end
 # Apply functions for diagonalising different shapes of transfer matrices
 # =======================================================================
 # Fig.25 of https://arxiv.org/pdf/2311.18785. Firstly appear in Chenfeng Bao's thesis, see http://hdl.handle.net/10012/14674.
+"""
+When the elementary modular parameter for TA, TB is `τ`,
+the transfer matrix has `τ_TM = (1 + τ) / 2 / (1 - τ)`.
+"""
 function MPO_action_2gates(TA::TensorMap, TB::TensorMap, x::TensorMap)
     @tensor fx[-1 -2 -3 -4; 5] := TB[-1 -2; 1 2] * x[1 2 3 4; 5] * TB[-3 -4; 3 4]
     @tensor ffx[-1 -2 -3 -4; 5] := TA[-3 -4; 2 3] * fx[1 2 3 4; 5] *
@@ -259,6 +282,10 @@ function MPO_action_2gates(TA::TensorMap, TB::TensorMap, x::TensorMap)
     return permute(ffx, ((2, 3, 4, 1), (5,)))
 end
 
+"""
+When the elementary modular parameter for TA, TB is `τ`,
+the transfer matrix has `τ_TM = τ / 4`.
+"""
 function MPO_action_1x4(TA::TensorMap, TB::TensorMap, x::TensorMap)
     @tensor TTTTx[-1 -2 -3 -4; -5] := x[1 2 3 4; -5] * TA[41 -1; 1 12] *
         TB[12 -2; 2 23] *
@@ -266,18 +293,25 @@ function MPO_action_1x4(TA::TensorMap, TB::TensorMap, x::TensorMap)
     return TTTTx
 end
 
+"""
+When the elementary modular parameter for TA, TB is `τ`,
+the transfer matrix has `τ_TM = (τ + 1) / 4`.
+"""
 function MPO_action_1x4_twist(TA::TensorMap, TB::TensorMap, x::TensorMap)
     TTTTx = MPO_action_1x4(TA, TB, x)
     return permute(TTTTx, ((2, 3, 4, 1), (5,)))
 end
 
+"""
+This renormalization will change the elementary
+modular parameter from τ to (τ - 1) / 2.
+"""
 function reduced_MPO(
         dl::TensorMap, ur::TensorMap, ul::TensorMap, dr::TensorMap,
         trunc::TruncationStrategy
     )
-    @plansor temp[-1 -2; -3 -4] := ur[-1; 1 4] *
-        ul[4; 3 -2] *
-        dr[-3; 2 1] * dl[2; -4 3]
+    @plansor temp[-1 -2; -3 -4] :=
+        ur[-1; 1 4] * ul[4; 3 -2] * dr[-3; 2 1] * dl[2; -4 3]
     D, U = SVD12(temp, trunc)
     @plansor translate[-1 -2; -3 -4] := U[-2; 1 -4] * D[-1 1; -3]
     return translate

From 03b35f07a0401adba87dbd40214a4b714fcf6d14 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Tue, 9 Jun 2026 20:47:34 +0800
Subject: [PATCH 08/14] Update loop-TNR test with anisotropic Ising model

---
 test/schemes/schemes.jl | 50 ++++++++++++++++++++++++-----------------
 1 file changed, 30 insertions(+), 20 deletions(-)

diff --git a/test/schemes/schemes.jl b/test/schemes/schemes.jl
index 70d4efed..fe841a11 100644
--- a/test/schemes/schemes.jl
+++ b/test/schemes/schemes.jl
@@ -166,10 +166,15 @@ end
     @test X2 ≈ 2.0 rtol = 1.0e-2
 end
 
-# LoopTNR
-@testset "LoopTNR - Ising Model - Dense Solver" begin
-    @info "LoopTNR ising free energy"
-    scheme = LoopTNR(T)
+# LoopTNR — anisotropic Ising
+const Jx_aniso, Jy_aniso = 1.0, 0.6
+const βc_aniso_test = ising_anisotropic_βc(Jx_aniso, Jy_aniso)
+const T_aniso = classical_ising(Z2Irrep, βc_aniso_test; Jx = Jx_aniso, Jy = Jy_aniso)
+const f_aniso_exact = f_onsager_anisotropic(βc_aniso_test, Jx_aniso, Jy_aniso)
+
+@testset "LoopTNR - Anisotropic Ising Model (Jx=1.0, Jy=0.6)" begin
+    @info "LoopTNR anisotropic ising free energy"
+    scheme = LoopTNR(T_aniso)
 
     loop_condition = LoopParameters(
         sweeping = maxiter(5) & convcrit(1.0e-9, (steps, cost) -> abs(cost[end])),
@@ -180,30 +185,35 @@ end
         scheme, truncrank(8), maxiter(25), loop_condition
     )
 
-    @test free_energy(data, ising_βc) ≈ f_onsager rtol = 1.0e-6
+    @test free_energy(data, βc_aniso_test) ≈ f_aniso_exact rtol = 1.0e-6
 
-    @info "LoopTNR ising CFT data"
-    scheme = LoopTNR(T)
+    @info "LoopTNR anisotropic ising CFT data"
+    scheme = LoopTNR(T_aniso)
     run!(scheme, truncrank(12), maxiter(10))
 
     for shape in [[1, 4, 1], [sqrt(2), 2 * sqrt(2), 0]]
-        cft = sort(CFTData(scheme; shape = shape).scaling_dimensions; by = abs)
-        d1, d2 = real(cft[Z2Irrep(1)][1]), real(cft[Z2Irrep(0)][2])
-        @info "Obtained lowest scaling dimensions:\n$(d1), $(d2)."
-        @test d1 ≈ ising_cft_exact[1] rtol = 5.0e-4
-        @test d2 ≈ ising_cft_exact[2] rtol = 5.0e-4
+        cft = CFTData(scheme; shape = shape)
+        d_σ = real(cft.scaling_dimensions[Z2Irrep(1)][1])
+        d_ε = real(cft.scaling_dimensions[Z2Irrep(0)][2])
+        @info "Shape $shape:  Δ(σ) = $d_σ,  Δ(ε) = $d_ε,  c = $(cft.central_charge)"
+        @test d_σ ≈ ising_cft_exact[1] rtol = 5.0e-4
+        @test d_ε ≈ ising_cft_exact[2] rtol = 5.0e-4
+        @test cft.central_charge ≈ 0.5 rtol = 5.0e-3
     end
 
     for shape in [[1, 8, 1], [4 / sqrt(10), 2 * sqrt(10), 2 / sqrt(10)]]
-        cft = sort(CFTData(scheme; shape = shape, trunc = truncrank(16), truncentanglement = trunctol(atol = 1.0e-10)).scaling_dimensions; by = real)
-        d1, d2 = real(cft[Z2Irrep(1)][1]), real(cft[Z2Irrep(0)][2])
-        @info "Obtained lowest scaling dimensions:\n$(d1), $(d2)."
-        @test d1 ≈ ising_cft_exact[1] rtol = 1.0e-3
-        @test d2 ≈ ising_cft_exact[2] rtol = 1.0e-3
+        cft = CFTData(scheme; shape = shape, trunc = truncrank(16),
+                      truncentanglement = trunctol(atol = 1.0e-10))
+        d_σ = real(cft.scaling_dimensions[Z2Irrep(1)][1])
+        d_ε = real(cft.scaling_dimensions[Z2Irrep(0)][2])
+        @info "Shape $shape:  Δ(σ) = $d_σ,  Δ(ε) = $d_ε,  c = $(cft.central_charge)"
+        @test d_σ ≈ ising_cft_exact[1] rtol = 2.0e-3
+        @test d_ε ≈ ising_cft_exact[2] rtol = 2.0e-3
+        @test cft.central_charge ≈ 0.5 rtol = 1.0e-2
     end
 
-    @info "LoopTNR ising ground state degeneracy"
-    T1 = classical_ising(ising_βc - 0.01)
+    @info "LoopTNR anisotropic ising ground state degeneracy"
+    T1 = classical_ising(βc_aniso_test - 0.01; Jx = Jx_aniso, Jy = Jy_aniso)
     scheme = LoopTNR(T1)
     run!(scheme, truncrank(12), maxiter(20))
     gsd = ground_state_degeneracy(scheme)
@@ -212,7 +222,7 @@ end
     @test X1 ≈ 1.0 rtol = 1.0e-2
     @test X2 ≈ 1.0 rtol = 1.0e-2
 
-    T2 = classical_ising(ising_βc + 0.01)
+    T2 = classical_ising(βc_aniso_test + 0.01; Jx = Jx_aniso, Jy = Jy_aniso)
     scheme = LoopTNR(T2)
     run!(scheme, truncrank(12), maxiter(20))
     gsd = ground_state_degeneracy(scheme)

From 50e5ceadec6bf63d4c0ef73b4e89182865d1e4ee Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Wed, 10 Jun 2026 09:51:22 +0800
Subject: [PATCH 09/14] Fix merge regressions

---
 src/utility/cft.jl      | 6 +++---
 test/schemes/schemes.jl | 6 ++++--
 2 files changed, 7 insertions(+), 5 deletions(-)

diff --git a/src/utility/cft.jl b/src/utility/cft.jl
index c88c8e04..4858f89f 100644
--- a/src/utility/cft.jl
+++ b/src/utility/cft.jl
@@ -1,5 +1,5 @@
 """
-    CFTData{E, I} where {E, I}
+    struct CFTData{E, K, V, A <: AbstractVector{E}}
 
 A struct to hold conformal data extracted from a TNR scheme.
 
@@ -11,10 +11,10 @@ A struct to hold conformal data extracted from a TNR scheme.
 # Fields
     - `central_charge::Union{E, Missing}`: The central charge of the CFT. Will be `nothing` if not calculated.
     - `modular_parameter::E`: The elementary modular parameter of a square spacetime patch of the CFT.
-    - `scaling_dimensions::TensorKit.SectorVector{E, I}`: The scaling dimensions of the CFT, organized in a `TensorKit.SectorVector` 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, 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
 
 """
-struct CFTData{E, I}
+struct CFTData{E, K, V, A <: AbstractVector{E}}
     "Central charge of the CFT. Will be `missing` if not calculated."
     central_charge::Union{E, Missing}
     "Elementary modular parameter for one tensor"
diff --git a/test/schemes/schemes.jl b/test/schemes/schemes.jl
index 67eaa6a9..2163de2e 100644
--- a/test/schemes/schemes.jl
+++ b/test/schemes/schemes.jl
@@ -202,8 +202,10 @@ const f_aniso_exact = f_onsager_anisotropic(βc_aniso_test, Jx_aniso, Jy_aniso)
     end
 
     for shape in [[1, 8, 1], [4 / sqrt(10), 2 * sqrt(10), 2 / sqrt(10)]]
-        cft = CFTData(scheme; shape = shape, trunc = truncrank(16),
-                      truncentanglement = trunctol(atol = 1.0e-10))
+        cft = CFTData(
+            scheme; shape = shape, trunc = truncrank(16),
+            truncentanglement = trunctol(atol = 1.0e-10)
+        )
         d_σ = real(cft.scaling_dimensions[Z2Irrep(1)][1])
         d_ε = real(cft.scaling_dimensions[Z2Irrep(0)][2])
         @info "Shape $shape:  Δ(σ) = $d_σ,  Δ(ε) = $d_ε,  c = $(cft.central_charge)"

From eeead5714827546040248e453df72a8c66257eae Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Wed, 10 Jun 2026 10:16:13 +0800
Subject: [PATCH 10/14] Fix [1, 1, 0] and test with TRG

---
 src/utility/cft.jl      | 15 ++++++++++-----
 test/schemes/schemes.jl | 39 +++++++++++++++++++++------------------
 2 files changed, 31 insertions(+), 23 deletions(-)

diff --git a/src/utility/cft.jl b/src/utility/cft.jl
index 4858f89f..f82825d9 100644
--- a/src/utility/cft.jl
+++ b/src/utility/cft.jl
@@ -42,7 +42,7 @@ end
 function CFTData(T::TensorMap{E, S, 2, 2}; shape = [sqrt(2), 2 * sqrt(2), 0], kwargs...) where {E, S}
     if shape == [1, 1, 0] # trivial implementation
         τ0 = elementary_modular_parameter(T, T)
-        Δs = _scaling_dimensions(T)
+        Δs = _scaling_dimensions(T, τ0)
         return CFTData(missing, τ0, StructuredVector(Δs, Dict([Trivial => collect(eachindex(Δs))])))
     else
         CFTData(T, T; shape, kwargs...)
@@ -70,9 +70,12 @@ function CFTData(
     end
 end
 
-# TODO: replace v with solved elementary modular parameter
-function _scaling_dimensions(T::TensorMap{E, S, 2, 2}; v = 1, unitcell = 1) where {E, S}
-    # stack unitcell copies of T and trace
+"""
+Construct the transfer matrix along vertical direction
+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
 
@@ -88,7 +91,9 @@ function _scaling_dimensions(T::TensorMap{E, S, 2, 2}; v = 1, unitcell = 1) wher
     data = filter(x -> real(x) > 0, data) # filtering out negative real values
     data = filter(x -> abs(x) > 1.0e-12, data) # filtering out small values
 
-    return unitcell * (1 / (2π * v)) .* log.(data[1] ./ data)
+    # modular parameter of the constructed transfer matrix
+    Imτ = imag(τ0) / unitcell
+    return 1 / (2π * Imτ) .* log.(data[1] ./ data)
 end
 
 """
diff --git a/test/schemes/schemes.jl b/test/schemes/schemes.jl
index 2163de2e..1893a544 100644
--- a/test/schemes/schemes.jl
+++ b/test/schemes/schemes.jl
@@ -11,31 +11,39 @@ T = classical_ising()
 T_3D = classical_ising_3D()
 const f_benchmark3D = ising_3D_free_energy_htse()
 
+const Jx_aniso, Jy_aniso = 1.0, 0.6
+const βc_aniso_test = ising_anisotropic_βc(Jx_aniso, Jy_aniso)
+const T_aniso = classical_ising(Z2Irrep, βc_aniso_test; Jx = Jx_aniso, Jy = Jy_aniso)
+const f_aniso_exact = f_onsager_anisotropic(βc_aniso_test, Jx_aniso, Jy_aniso)
+
 function cft_finalize!(scheme)
     finalize!(scheme)
     return CFTData(scheme)
 end
 
 # TRG
-@testset "TRG - Ising Model" begin
-    @info "TRG ising free energy"
-    scheme = TRG(T)
+@testset "TRG - Anisotropic Ising Model (Jx=1.0, Jy=0.6)" begin
+    @info "TRG anisotropic ising free energy"
+    scheme = TRG(T_aniso)
     data = run!(scheme, truncrank(24), maxiter(25))
 
-    @test free_energy(data, ising_βc) ≈ f_onsager rtol = 2.0e-6
+    @test free_energy(data, βc_aniso_test) ≈ f_aniso_exact rtol = 2.0e-6
 
-    @info "TRG ising CFT data"
-    scheme = TRG(T)
+    @info "TRG anisotropic ising CFT data — shape [1, 1, 0]"
+    scheme = TRG(T_aniso)
     run!(scheme, truncrank(24), maxiter(10))
 
-    cft = sort(CFTData(scheme; shape = [1, 1, 0]).scaling_dimensions[2:end]; by = abs) .|> real
+    cft = CFTData(scheme; shape = [1, 1, 0])
+    sd_all = real(cft.scaling_dimensions[Trivial])
+    cft_sorted = sort(sd_all[2:end]; by = abs)
 
-    @test cft[1] ≈ ising_cft_exact[1] rtol = 2.0e-4
-    @test cft[2] ≈ ising_cft_exact[2] rtol = 1.0e-2
+    @test cft_sorted[1] ≈ ising_cft_exact[1] rtol = 2.0e-4
+    @test cft_sorted[2] ≈ ising_cft_exact[2] rtol = 2.0e-2
 
-    @info "TRG ising ground state degeneracy"
+    @info "Obtained scaling dimensions: Δ₁ = $(cft_sorted[1]), Δ₂ = $(cft_sorted[2])"
 
-    T1 = classical_ising(ising_βc - 0.01)
+    @info "TRG anisotropic ising ground state degeneracy"
+    T1 = classical_ising(βc_aniso_test - 0.01; Jx = Jx_aniso, Jy = Jy_aniso)
     scheme = TRG(T1)
     run!(scheme, truncrank(16), maxiter(20))
     gsd = ground_state_degeneracy(scheme)
@@ -44,7 +52,7 @@ end
     @test X1 ≈ 1.0 rtol = 1.0e-2
     @test X2 ≈ 1.0 rtol = 1.0e-2
 
-    T2 = classical_ising(ising_βc + 0.01)
+    T2 = classical_ising(βc_aniso_test + 0.01; Jx = Jx_aniso, Jy = Jy_aniso)
     scheme = TRG(T2)
     run!(scheme, truncrank(16), maxiter(20))
     gsd = ground_state_degeneracy(scheme)
@@ -52,7 +60,6 @@ end
     @test gsd ≈ 2 rtol = 1.0e-2
     @test X1 ≈ 2.0 rtol = 1.0e-2
     @test X2 ≈ 2.0 rtol = 1.0e-2
-
 end
 
 # BTRG
@@ -166,11 +173,7 @@ end
     @test X2 ≈ 2.0 rtol = 1.0e-2
 end
 
-# LoopTNR — anisotropic Ising
-const Jx_aniso, Jy_aniso = 1.0, 0.6
-const βc_aniso_test = ising_anisotropic_βc(Jx_aniso, Jy_aniso)
-const T_aniso = classical_ising(Z2Irrep, βc_aniso_test; Jx = Jx_aniso, Jy = Jy_aniso)
-const f_aniso_exact = f_onsager_anisotropic(βc_aniso_test, Jx_aniso, Jy_aniso)
+# LoopTNR
 
 @testset "LoopTNR - Anisotropic Ising Model (Jx=1.0, Jy=0.6)" begin
     @info "LoopTNR anisotropic ising free energy"

From b2a843fc134d5eca4cdcf35cbcac1af90a1d766e Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Wed, 10 Jun 2026 11:23:50 +0800
Subject: [PATCH 11/14] Docstring improvement

---
 src/utility/cft.jl | 18 ++++++++----------
 1 file changed, 8 insertions(+), 10 deletions(-)

diff --git a/src/utility/cft.jl b/src/utility/cft.jl
index f82825d9..1bd84a21 100644
--- a/src/utility/cft.jl
+++ b/src/utility/cft.jl
@@ -10,7 +10,7 @@ A struct to hold conformal data extracted from a TNR scheme.
 
 # Fields
     - `central_charge::Union{E, Missing}`: The central charge of the CFT. Will be `nothing` if not calculated.
-    - `modular_parameter::E`: The elementary modular parameter of a square spacetime patch of the CFT.
+    - `modular_parameter::E`: The elementary modular parameter of a square patch of the lattice network.
     - `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
 
 """
@@ -120,18 +120,16 @@ end
 
 # The case with spin is based on https://arxiv.org/pdf/1512.03846 and some private communications with Yingjie Wei and Atsushi Ueda
 """
-    spec(TA::TensorMap, TB::TensorMap, shape::Array; Nh = 25)
-
 Internal function to construct transfer matrices and extract conformal data.
 
-# Parameters
-
-TA, TB: Tensors used to construct the transfer matrix. 
-    They can be different from the original tensor in the network.
-shape: A triplet `[h, L, x]` describing the geometry of the transfer matrix.
-τ0: The modular parameter of a square patch of the original network.
+# Arguments
+- `TA, TB`: Rank-4 network tensors (may be identical for 1-site unit cells).
+- `shape`:  A triplet `[h, L, x]` — height, circumference, and horizontal shift
+  of the tube geometry, in units of the original tensor patch.
+- `τ0`:     Elementary modular parameter of one tensor.
+- `Nh`:     Number of eigenvalues to solve for per sector (default 25).
 """
-function spec(TA::TensorMap, TB::TensorMap, shape::Array, τ0::Number; Nh = 25)
+function spec(TA::TensorMap, TB::TensorMap, shape::Vector{<:Number}, τ0::Number; Nh = 25)
     I = sectortype(TA)
     𝔽 = field(TA)
     if BraidingStyle(I) != Bosonic()

From db393a3552ecf4e84bd7749b8cc30d9c2a89fd54 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Wed, 10 Jun 2026 21:12:00 +0800
Subject: [PATCH 12/14] Better way to solve tau with only the leading
 eigenvalues

---
 Project.toml       |  2 +
 src/TNRKit.jl      |  1 +
 src/utility/cft.jl | 95 ++++++++++++++++++++++++++++++++++++----------
 3 files changed, 77 insertions(+), 21 deletions(-)

diff --git a/Project.toml b/Project.toml
index b0ac7086..2a380006 100644
--- a/Project.toml
+++ b/Project.toml
@@ -11,6 +11,7 @@ KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LoggingExtras = "e6f89c97-d47a-5376-807f-9c37f3926c36"
 MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4"
+NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
@@ -28,6 +29,7 @@ FastGaussQuadrature = "1.1.0"
 KrylovKit = "0.10"
 LoggingExtras = "~1.0"
 MatrixAlgebraKit = "0.6.1"
+NonlinearSolve = "4.12.0"
 OptimKit = "0.4.0"
 QuadGK = "2.11.2"
 Roots = "3.0.0"
diff --git a/src/TNRKit.jl b/src/TNRKit.jl
index 3c6243d8..53a0ebc6 100644
--- a/src/TNRKit.jl
+++ b/src/TNRKit.jl
@@ -11,6 +11,7 @@ using SpecialFunctions
 using FastGaussQuadrature
 using QuadGK
 using Roots
+using NonlinearSolve
 using Base.Threads
 using Combinatorics: permutations
 
diff --git a/src/utility/cft.jl b/src/utility/cft.jl
index 1bd84a21..1c16da6c 100644
--- a/src/utility/cft.jl
+++ b/src/utility/cft.jl
@@ -291,23 +291,24 @@ Extract the elementary modular parameter of one tensor from 2x2 transfer matrice
 function elementary_modular_parameter(
         TA::TensorMap{E, S, 2, 2}, TB::TensorMap{E, S, 2, 2}
     ) where {E, S}
-    # vertical (north to south)
-    evs_v = _eigsolve_2x2_NtoS(TA, TB)
-    # horizontal (east to west)
-    evs_h = _eigsolve_2x2_EtoW(TA, TB)
-    # diagonal (northeast to southwest)
-    evs_a = _eigsolve_2x2_NEtoSW(TA, TB)
-    # diagonal (northwest to southeast)
-    evs_b = _eigsolve_2x2_NWtoSE(TA, TB)
-    # norm²
-    v = sqrt(
-        log(abs(evs_v[2] / evs_v[1])) /
-            log(abs(evs_h[2] / evs_h[1]))
-    )
-    # phase
-    r = log(abs(evs_a[2] / evs_a[1])) /
-        log(abs(evs_b[2] / evs_b[1]))
-    θ = acos((v^2 + 1) / (2 * v) * (r - 1) / (r + 1))
+    # TODO: fermions
+    if BraidingStyle(sectortype(TA)) != Bosonic
+        return complex(0.0, 1.0)
+    end
+    # leading eigenvalue of each transfer matrix
+    # corresponding to the Δ = s = 0 primary field
+    λv = _eigsolve_2x2_NtoS(TA, TB)
+    λh = _eigsolve_2x2_EtoW(TA, TB)
+    λa = _eigsolve_2x2_NEtoSW(TA, TB)
+    λb = _eigsolve_2x2_NWtoSE(TA, TB)
+    # edge case: c = 0
+    if all(isapprox.(λv, (λh, λa, λb); rtol = 1.0e-6))
+        return complex(0.0, 1.0)
+    end
+    # when c ≠ 0
+    a1, a2, a3 = log(λh / λv), log(λa / λv), log(λb / λv)
+    # c here is actually π c / 6
+    c, v, θ = solve_cvtheta(a1, a2, a3)
     return v * cis(θ)
 end
 
@@ -318,11 +319,12 @@ function _eigsolve_2x2_NtoS(TA, TB)
         return fx
     end
     x0 = ones(domain(TA, 1) ⊗ domain(TB, 1))
-    spec0, _, info = eigsolve(f0, x0, 2, :LM; verbosity = 0)
+    spec0, _, info = eigsolve(f0, x0, 1, :LM; verbosity = 0)
     if info.converged == 0
         @warn "The area term eigensolver did not converge."
     end
-    return spec0
+    # TODO: for fermions
+    return first(spec0)
 end
 
 function _eigsolve_2x2_EtoW(TA, TB)
@@ -350,11 +352,12 @@ function _eigsolve_2x2_NEtoSW(TA, TB)
         return fx
     end
     x0 = ones(domain(TA, 1) ⊗ domain(TB, 1) ⊗ domain(TB, 2) ⊗ domain(TA, 2))
-    spec0, _, info = eigsolve(f0, x0, 2, :LM; verbosity = 0)
+    spec0, _, info = eigsolve(f0, x0, 1, :LM; verbosity = 0)
     if info.converged == 0
         @warn "The area term eigensolver did not converge."
     end
-    return spec0
+    # TODO: for fermions
+    return first(spec0)
 end
 
 function _eigsolve_2x2_NWtoSE(TA, TB)
@@ -362,3 +365,53 @@ function _eigsolve_2x2_NWtoSE(TA, TB)
     TB′ = permute(TB, ((2, 4), (1, 3)))
     return _eigsolve_2x2_NEtoSW(TB′, TA′)
 end
+
+# Logistic sigmoid and its inverse (logit)
+sigmoid(x) = 1 / (1 + exp(-x))
+logit(p) = log(p / (1 - p))
+
+"""
+    solve_cvtheta(a1, a2, a3; c0=0.5, v0=1.0, θ0=π/2)
+
+Solve for positive (c, v) and θ ∈ (0, π) from the three equations:
+
+    a1 = (1/v - v) * c * sin(θ)
+    a2 = (v / (1 + v² - 2v cos θ) - v) * c * sin(θ)
+    a3 = (v / (1 + v² + 2v cos θ) - v) * c * sin(θ)
+
+Returns `(c, v, θ)` as a named tuple.
+"""
+function solve_cvtheta(a1, a2, a3; c0 = 0.5, v0 = 1.0, θ0 = π / 2)
+    # Work in unconstrained coords to keep variables in their natural domain.
+    #   c = exp(xc)        →  c > 0
+    #   v = exp(xv)        →  v > 0
+    #   θ = π * σ(xθ)      →  θ ∈ (0, π)
+    function f!(du, u, p)
+        xc, xv, xθ = u
+        c = exp(xc)
+        v = exp(xv)
+        θ = π * sigmoid(xθ)
+
+        sinθ = sin(θ)
+        cosθ = cos(θ)
+
+        du[1] = (1 / v - v) * c * sinθ - a1
+        du[2] = (v / (1 + v^2 - 2v * cosθ) - v) * c * sinθ - a2
+        return du[3] = (v / (1 + v^2 + 2v * cosθ) - v) * c * sinθ - a3
+    end
+
+    # Initial guess in unconstrained space
+    u0 = [log(c0), log(v0), logit(θ0 / π)]
+    prob = NonlinearProblem(f!, u0)
+    sol = solve(prob, NewtonRaphson(; autodiff = AutoForwardDiff()))
+
+    if !SciMLBase.successful_retcode(sol)
+        @warn "Solver did not converge" retcode = sol.retcode residuals = sol.residuals
+    end
+
+    xc, xv, xθ = sol.u
+    c = exp(xc)
+    v = exp(xv)
+    θ = π * sigmoid(xθ)
+    return c, v, θ
+end

From c7e9c88b3adbae7bd1a340864e5680016c8d9b8d Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Wed, 10 Jun 2026 21:21:04 +0800
Subject: [PATCH 13/14] Fix docstring

---
 src/utility/cft.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/utility/cft.jl b/src/utility/cft.jl
index 1c16da6c..175620b4 100644
--- a/src/utility/cft.jl
+++ b/src/utility/cft.jl
@@ -379,7 +379,7 @@ Solve for positive (c, v) and θ ∈ (0, π) from the three equations:
     a2 = (v / (1 + v² - 2v cos θ) - v) * c * sin(θ)
     a3 = (v / (1 + v² + 2v cos θ) - v) * c * sin(θ)
 
-Returns `(c, v, θ)` as a named tuple.
+Returns `(c, v, θ)`.
 """
 function solve_cvtheta(a1, a2, a3; c0 = 0.5, v0 = 1.0, θ0 = π / 2)
     # Work in unconstrained coords to keep variables in their natural domain.

From 957fe145b2cc757c85630346307d51748f2dbed8 Mon Sep 17 00:00:00 2001
From: Yue Zhengyuan <yuezy1997@icloud.com>
Date: Wed, 10 Jun 2026 22:46:13 +0800
Subject: [PATCH 14/14] Fix fermion check and c = 0 detection

---
 src/utility/cft.jl | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/utility/cft.jl b/src/utility/cft.jl
index 175620b4..cab36eb9 100644
--- a/src/utility/cft.jl
+++ b/src/utility/cft.jl
@@ -292,15 +292,15 @@ function elementary_modular_parameter(
         TA::TensorMap{E, S, 2, 2}, TB::TensorMap{E, S, 2, 2}
     ) where {E, S}
     # TODO: fermions
-    if BraidingStyle(sectortype(TA)) != Bosonic
+    if BraidingStyle(sectortype(TA)) != Bosonic()
         return complex(0.0, 1.0)
     end
     # leading eigenvalue of each transfer matrix
-    # corresponding to the Δ = s = 0 primary field
-    λv = _eigsolve_2x2_NtoS(TA, TB)
-    λh = _eigsolve_2x2_EtoW(TA, TB)
-    λa = _eigsolve_2x2_NEtoSW(TA, TB)
-    λb = _eigsolve_2x2_NWtoSE(TA, TB)
+    # corresponding to the Δ = s = 0 identity field
+    λv = real(_eigsolve_2x2_NtoS(TA, TB))
+    λh = real(_eigsolve_2x2_EtoW(TA, TB))
+    λa = real(_eigsolve_2x2_NEtoSW(TA, TB))
+    λb = real(_eigsolve_2x2_NWtoSE(TA, TB))
     # edge case: c = 0
     if all(isapprox.(λv, (λh, λa, λb); rtol = 1.0e-6))
         return complex(0.0, 1.0)