Fix the bug of modular condition of S matrix

src/sectors.jl

@@ -691,15 +691,15 @@ hopflink(a::I, b::I) where {I <: Sector} = sum(dim(c) * tr(Rsymbol(a, b, c) * Rs
 """
      Smatrix(::Type{I}) where {I <: Sector}
 
-Return the S-matrix of the sector type `I`, which is a matrix containing the hopflinks of all pairs of sectors of type `I`.
+Return the S-matrix of the sector type `I`, which is a matrix containing the hopflinks of all pairs of sectors of type `I`, with the second sector being taken dual.
 The S-matrix is not normalized by the total quantum dimension here.
 """
 function Smatrix(::Type{I}) where {I <: Sector}
     Base.IteratorSize(values(I)) isa Base.IsInfinite &&
         throw(ArgumentError("Only defined for sectors with a finite number of simple objects"))
     vals = values(I)
     l = length(vals)
-    return reshape([hopflink(a, b) for a in vals, b in vals], (l, l))
+    return reshape([hopflink(a, dual(b)) for a in vals, b in vals], (l, l))
 end
 
 """
@@ -715,13 +715,14 @@ end
 """
     topological_central_charge(::Type{I}) where {I <: Sector}
 
-Return the topological central charge c of the modular sector type `I`, where c is determined mod 8.
+Return the topological central charge c of the braided sector type `I`, where c is determined mod 8.
 We choose convention by restrict the returning value as rational numbers in (-4, 4].
 """
 function topological_central_charge(::Type{I}) where {I <: Sector}
-    ξ = sum(dim(a)^2 * twist(a) for a in values(I)) / dim(I)
-    @assert isapprox(abs(ξ), 1) "Sector $I is not modular"
-    c_float = angle(ξ) * 8 / (2π)
+    gauss_sum = sum(dim(a)^2 * twist(a) for a in values(I))
+    Θ = abs(gauss_sum)
+    @assert Θ > sqrt(eps(float(Θ))) "Topological central charge is not defined for sector type $I" # For non-modular categories, central charge is also meaningful. See https://arxiv.org/pdf/1602.05946. For super modular category, Gauss sum vanishes, and its central charge needs to be defined in another manner: https://arxiv.org/pdf/1603.09294.
+    c_float = angle(gauss_sum) * 8 / (2π)
 
     isapprox(c_float, -4) && return 4 // 1
 

test/runtests.jl

@@ -266,15 +266,22 @@ end
 end
 
 @testset "Ismodular" begin
-    @test !ismodular(Z2Irrep)
-    @test !ismodular(Z3Irrep)
-    @test !ismodular(FermionParity)
-    @test !ismodular(A4Irrep)
-    @test !ismodular(IsingAnyon ⊠ Z2Irrep)
-    @test ismodular(IsingAnyon)
-    @test ismodular(FibonacciAnyon)
-    @test ismodular(TimeReversed{IsingAnyon})
-    @test ismodular(IsingAnyon ⊠ TimeReversed{IsingAnyon})
+    Tannakian_list = [Z2Irrep, Z3Irrep, A4Irrep, D3Irrep, D4Irrep, Z2Irrep ⊠ D4Irrep, D3Irrep ⊠ A4Irrep]
+    Super_Tannakian_list = [FermionParity, FermionParity ⊠ A4Irrep, FermionParity ⊠ Z3Irrep, D4Irrep ⊠ FermionParity]
+    UMTC_list = [
+        IsingAnyon, FibonacciAnyon, TimeReversed{IsingAnyon}, TimeReversed{FibonacciAnyon},
+        FibonacciAnyon ⊠ FibonacciAnyon, FibonacciAnyon ⊠ IsingAnyon, IsingAnyon ⊠ TimeReversed{IsingAnyon},
+        TimeReversed{FibonacciAnyon} ⊠ IsingAnyon, IsingAnyon ⊠ IsingAnyon ⊠ IsingAnyon,
+        IsingAnyon ⊠ FibonacciAnyon ⊠ IsingAnyon ⊠ TimeReversed{IsingAnyon} ⊠ TimeReversed{FibonacciAnyon},
+    ]
+    UMTC_over_RepG_list = [Z2Irrep ⊠ IsingAnyon, Z3Irrep ⊠ FibonacciAnyon, D3Irrep ⊠ TimeReversed{IsingAnyon}, A4Irrep ⊠ FibonacciAnyon ⊠ TimeReversed{IsingAnyon}]
+    UMTC_over_sRepG_list = [Z2Irrep ⊠ FermionParity ⊠ IsingAnyon, FermionParity ⊠ Z3Irrep ⊠ FibonacciAnyon, D3Irrep ⊠ TimeReversed{IsingAnyon} ⊠ FermionParity, A4Irrep ⊠ FibonacciAnyon ⊠ FermionParity ⊠ TimeReversed{IsingAnyon}]
+    for sect in [Tannakian_list..., Super_Tannakian_list..., UMTC_over_RepG_list..., UMTC_over_sRepG_list...]
+        @test !ismodular(sect)
+    end
+    for sect in UMTC_list
+        @test ismodular(sect)
+    end
 end
 
 @testset "Total quantum dimension" begin
@@ -290,6 +297,10 @@ end
 end
 
 @testset "Topological central charge" begin
+    @test topological_central_charge(Z2Irrep) == 0 // 1
+    @test topological_central_charge(A4Irrep) == 0 // 1
+    @test topological_central_charge(Z2Irrep ⊠ IsingAnyon) == 1 // 2
+    @test topological_central_charge(FibonacciAnyon ⊠ D4Irrep) == - 14 // 5
     @test topological_central_charge(IsingAnyon) == 1 // 2
     @test topological_central_charge(TimeReversed{IsingAnyon}) == - 1 // 2
     @test topological_central_charge(IsingAnyon ⊠ TimeReversed{IsingAnyon}) == 0 // 1