Potts edits (#29)

* use new MPSKit function add_util_mpoleg for mpoham

* generalise Potts operators to general Q

* add Potts tests for Q=4 and Q=5

* Revert "use new MPSKit function add_util_mpoleg for mpoham"

This reverts commit ba4a409.

* generalise Potts operators to general Q

* add Potts tests for Q=4 and Q=5

* remove U_matrix, V_matrix, add weyl_heisenberg_matrices, potts_Z, potts_X

* fix typo in weyl_heisenberg_matrices

* fix kwargs for potts_Z, potts_X

* fix typo in weyl_heisenberg_matrices: the sequel

* update potts_exchange, potts_field for Trivial symmetry

* update potts_exchange, potts_field for Z_Q symmetry

* clean up theoretical energy function

* redefine potts_exchange to be the ZZ' term of the hamiltonian, i.e., remove sum

* redefine potts_field to be fully dependent on potts_X, added potts_X with Z_Q symmetry

* add TODO for potts_Z with Z_Q symmetry (currently uses potts_X)

* update Potts hamiltonian to perform extra sum

* construct potts_X with Z_Q symmetry via blocks instead of array data input

* add potts_Z with Z_Q symmetry, update potts_Z to use this symmetric potts_Z

* export weyl_heisenberg_matrices, potts_X, potts_Z, potts_exchange, potts_field

* fix definition of W matrix in weyl_heisenberg_matrices

* rename potts_exchange to potts_ZZ, construct potts_ZZ with Z_Q symmetry via block structure

* remove old potts_ZZ with Z_Q symmetry, remove potts_Z with Z_Q symmetry since ill-defined

* update quantum_potts to use potts_ZZ

* add potts_ZZ to export, remove potts_exchange

* add tests for non-symmetric and Z3 symmetric 3-state Potts

* update tests for Z_Q symmetric potts operators

* edit description for potts_field

* add documentation for potts operators

* julia formatter

* use add_util_mpoleg in LocalOperator

* Revert "use add_util_mpoleg in LocalOperator"

This reverts commit 3a9adf2.

Author: borisdevos

docs/src/man/operators.md

@@ -52,6 +52,18 @@ For convenience, the Pauli matrices can also be recovered as ``σⁱ = 2 Sⁱ``.
 σσ
 ```
 
+## Q-state Potts operators
+
+The Q-state Potts operators `potts_X` and `potts_Z` are defined to fulfill the braiding relation ``ZX = \omega XZ`` with ``\omega = e^{2\pi i/Q}``. 
+
+Supported values of `symmetry` for the `X` operator are `Trivial` and `ZNIrrep{Q}`, while for the `Z` operator only `Trivial` is supported.
+
+```@docs
+potts_X
+potts_Z
+potts_ZZ
+```
+
 ## Bosonic operators
 
 The bosonic creation and annihilation operators `a_plus` ($$a^+$$) and `a_min` ($$a^-$$) are defined such that the following holds:

src/MPSKitModels.jl

@@ -19,12 +19,13 @@ export SnakePattern, frontandback_pattern, backandforth_pattern
 export LocalOperator, SumOfLocalOperators
 export @mpoham
 
-export spinmatrices, nonsym_spintensors, nonsym_bosonictensors
+export spinmatrices, nonsym_spintensors, nonsym_bosonictensors, weyl_heisenberg_matrices
 
 export S_x, S_y, S_z, S_plus, S_min
 export S_xx, S_yy, S_zz, S_plusmin, S_minplus, S_exchange
 export Sˣ, Sʸ, Sᶻ, S⁺, S⁻, Sˣˣ, Sʸʸ, Sᶻᶻ, S⁺⁻, S⁻⁺, SS
 export σˣ, σʸ, σᶻ, σ⁺, σ⁻, σˣˣ, σʸʸ, σᶻᶻ, σ⁺⁻, σ⁻⁺, σσ
+export potts_X, potts_Z, potts_field, potts_ZZ
 
 export a_plus, a_min, a_plusmin, a_minplus, a_number
 export a⁺, a⁻

src/models/hamiltonians.jl

@@ -264,11 +264,11 @@ function quantum_potts(elt::Type{<:Number}=ComplexF64,
                        symmetry::Type{<:Sector}=Trivial,
                        lattice::AbstractLattice=InfiniteChain(1);
                        q=3, J=1.0, g=1.0)
-    return @mpoham sum(nearest_neighbours(lattice)) do (i, j)
-        return -J * potts_exchange(elt, symmetry; q){i,j}
-    end - sum(vertices(lattice)) do i
-          return g * potts_field(elt, symmetry; q){i}
-          end
+    return @mpoham sum(sum(nearest_neighbours(lattice)) do (i, j)
+                           return -J * (potts_ZZ(elt, symmetry; q)^k){i,j}
+                       end - sum(vertices(lattice)) do i
+                             return g * (potts_field(elt, symmetry; q)^k){i}
+                             end for k in 1:(q - 1))
 end
 
 #===========================================================================================

src/operators/spinoperators.jl

@@ -446,62 +446,106 @@ const SS = S_exchange
 σσ(args...; kwargs...) = 4 * S_exchange(args...; kwargs...)
 
 """
-    potts_exchange([eltype::Type{<:Number}], [symmetry::Type{<:Sector}]; q=3)
+    potts_ZZ([eltype::Type{<:Number}], [symmetry::Type{<:Sector}]; q=3)
 
-The Potts exchange operator ``Z ⊗ Z' + Z' ⊗ Z``, where ``Z^q = 1``.
+The Potts operator ``Z ⊗ Z'``, where ``Z^q = 1``.
 """
-function potts_exchange end
-potts_exchange(; kwargs...) = potts_exchange(ComplexF64, Trivial; kwargs...)
-potts_exchange(elt::Type{<:Number}; kwargs...) = potts_exchange(elt, Trivial; kwargs...)
-function potts_exchange(symmetry::Type{<:Sector}; kwargs...)
-    return potts_exchange(ComplexF64, symmetry; kwargs...)
+function potts_ZZ end
+potts_ZZ(; kwargs...) = potts_ZZ(ComplexF64, Trivial; kwargs...)
+potts_ZZ(elt::Type{<:Number}; kwargs...) = potts_ZZ(elt, Trivial; kwargs...)
+function potts_ZZ(symmetry::Type{<:Sector}; kwargs...)
+    return potts_ZZ(ComplexF64, symmetry; kwargs...)
 end
 
-function potts_exchange(elt::Type{<:Number}, ::Type{Trivial}; q=3)
-    pspace = ComplexSpace(q)
-    Z = TensorMap(zeros, elt, pspace ← pspace)
-    for i in 1:q
-        Z[i, i] = cis(2π * (i - 1) / q)
-    end
-    return Z ⊗ Z' + Z' ⊗ Z
+function potts_ZZ(elt::Type{<:Number}, ::Type{Trivial}; q=3)
+    Z = potts_Z(elt, Trivial; q=q)
+    return Z ⊗ Z'
 end
-function potts_exchange(elt::Type{<:Number}, ::Type{ZNIrrep{Q}}; q=Q) where {Q}
+
+function potts_ZZ(elt::Type{<:Number}, ::Type{ZNIrrep{Q}}; q=Q) where {Q}
     @assert q == Q "q must match the irrep charge"
-    pspace = Vect[ZNIrrep{q}](i => 1 for i in 0:(q - 1))
-    aspace = Vect[ZNIrrep{q}](1 => 1, -1 => 1)
-    Z_left = TensorMap(ones, elt, pspace ← pspace ⊗ aspace)
-    Z_right = TensorMap(ones, elt, aspace ⊗ pspace ← pspace)
-    return contract_twosite(Z_left, Z_right)
+    pspace = Vect[ZNIrrep{Q}](i => 1 for i in 0:(Q - 1))
+    ZZ = TensorMap(zeros, elt, pspace ⊗ pspace ← pspace ⊗ pspace)
+    for charge in 0:(Q - 1)
+        for i in 1:Q
+            blocks(ZZ)[ZNIrrep{Q}(charge)][i, mod1(i + 1, Q)] = one(elt)
+        end
+    end
+    return ZZ
 end
 
 """
     potts_field([eltype::Type{<:Number}], [symmetry::Type{<:Sector}]; q=3) 
 
-The Potts field operator ``X + X'``, where ``X^q = 1``.
+The Potts field operator ``X``, an alias for ``potts_X``.
 """
-function potts_field end
-potts_field(; kwargs...) = potts_field(ComplexF64, Trivial; kwargs...)
-potts_field(elt::Type{<:Number}; kwargs...) = potts_field(elt, Trivial; kwargs...)
-function potts_field(symmetry::Type{<:Sector}; kwargs...)
-    return potts_field(ComplexF64, symmetry; kwargs...)
-end
+potts_field(args...; kwargs...) = potts_X(args...; kwargs...)
 
-function potts_field(elt::Type{<:Number}, ::Type{Trivial}; q=3)
-    pspace = ComplexSpace(q)
-    X = TensorMap(zeros, elt, pspace ← pspace)
-    for i in 1:q
-        X[mod1(i - 1, q), i] = one(elt)
+# Generalisations of Pauli matrices
+
+"""
+    weyl_heisenberg_matrices(dimension [, eltype])
+
+the Weyl-Heisenberg matrices according to [Wikipedia](https://en.wikipedia.org/wiki/Generalizations_of_Pauli_matrices#Sylvester's_generalized_Pauli_matrices_(non-Hermitian)).
+"""
+
+function weyl_heisenberg_matrices(Q::Int, elt=ComplexF64)
+    U = zeros(elt, Q, Q) # clock matrix
+    V = zeros(elt, Q, Q) # shift matrix
+    W = zeros(elt, Q, Q) # DFT
+    ω = cis(2 * pi / Q)
+
+    for row in 1:Q
+        U[row, row] = ω^(row - 1)
+        V[row, mod1(row - 1, Q)] = one(elt)
+        for col in 1:Q
+            W[row, col] = ω^((row - 1) * (col - 1))
+        end
     end
-    return X + X'
+    return U, V, W / sqrt(Q)
+end
+
+"""
+    potts_Z([eltype::Type{<:Number}], [symmetry::Type{<:Sector}]; Q=3)
+
+The Potts Z operator, also known as the clock operator, where Z^q=1.
+"""
+
+function potts_Z end
+potts_Z(; kwargs...) = potts_Z(ComplexF64, Trivial; kwargs...)
+potts_Z(elt::Type{<:Complex}; kwargs...) = potts_Z(elt, Trivial; kwargs...)
+potts_Z(symm::Type{<:Sector}; kwargs...) = potts_Z(ComplexF64, symm; kwargs...)
+
+function potts_Z(elt::Type{<:Number}, ::Type{Trivial}; q=3)
+    U, _, _ = weyl_heisenberg_matrices(q, elt)
+    Z = TensorMap(U, ComplexSpace(q) ← ComplexSpace(q))
+    return Z
 end
-# TODO: generalize to arbitrary q
-function potts_field(elt::Type{<:Number}, ::Type{ZNIrrep{Q}}; q=Q) where {Q}
+
+"""
+    potts_X([eltype::Type{<:Number}], [symmetry::Type{<:Sector}]; Q=3)
+
+The Potts X operator, also known as the shift operator, where X^q=1.
+"""
+
+function potts_X end
+potts_X(; kwargs...) = potts_X(ComplexF64, Trivial; kwargs...)
+potts_X(elt::Type{<:Complex}; kwargs...) = potts_X(elt, Trivial; kwargs...)
+potts_X(symm::Type{<:Sector}; kwargs...) = potts_X(ComplexF64, symm; kwargs...)
+
+function potts_X(elt::Type{<:Number}, ::Type{Trivial}; q=3)
+    _, V, _ = weyl_heisenberg_matrices(q, elt)
+    X = TensorMap(V, ComplexSpace(q) ← ComplexSpace(q))
+    return X
+end
+
+function potts_X(elt::Type{<:Number}, ::Type{ZNIrrep{Q}}; q=Q) where {Q}
     @assert q == Q "q must match the irrep charge"
-    @assert q == 3 "only q = 3 is implemented"
-    pspace = Vect[ZNIrrep{q}](i => 1 for i in 0:(q - 1))
+    pspace = Vect[ZNIrrep{Q}](i => 1 for i in 0:(Q - 1))
     X = TensorMap(zeros, elt, pspace ← pspace)
-    for (c, b) in blocks(X)
-        b .= isone(c) ? 2 : -1
+    ω = cis(2 * pi / Q)
+    for i in 1:Q
+        blocks(X)[ZNIrrep{Q}(i)] .= ω^i
     end
     return X
 end

test/potts.jl

@@ -2,20 +2,65 @@ using MPSKit
 using TensorKit
 
 alg = VUMPS(; maxiter=25, verbosity=0)
-E₀ = -(4 / 3 + 2sqrt(3) / π)
 
-@testset "Trivial" begin
+# https://iopscience.iop.org/article/10.1088/0305-4470/14/11/020/meta
+function E₀(Q::Int, maxiter::Int=1000)
+    Q == 3 && return -(4 / 3 + 2sqrt(3) / π)
+    Q == 4 && return 2 - 8 * log(2)
+    summation = sum((-1)^n / (sqrt(Q) / 2 - cosh((2 * n + 1) * acosh(sqrt(Q) / 2)))
+                    for n in 1:maxiter)
+    limit = 2 - Q - sqrt(Q) * (Q - 4) * summation
+    return limit
+end
+
+# Q = 3
+@testset "Q=3 Trivial" begin
     H = quantum_potts(; q=3)
     ψ = InfiniteMPS(3, 32)
     @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1e-10
     ψ, envs, δ = find_groundstate(ψ, H, alg)
-    @test E₀ ≈ expectation_value(ψ, H, envs) atol = 1e-2
+    @test E₀(3) ≈ expectation_value(ψ, H, envs) atol = 1e-2
 end
 
 @testset "Z3Irrep" begin
     H = quantum_potts(Z3Irrep; q=3)
     ψ = InfiniteMPS(Rep[ℤ₃](i => 1 for i in 0:2), Rep[ℤ₃](i => 12 for i in 0:2))
     @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1e-10
     ψ, envs, δ = find_groundstate(ψ, H, alg)
-    @test E₀ ≈ expectation_value(ψ, H, envs) atol = 1e-2
+    @test E₀(3) ≈ expectation_value(ψ, H, envs) atol = 1e-2
+end
+
+# Q = 4
+@testset "Q=4 Trivial" begin
+    H = quantum_potts(; q=4)
+    ψ = InfiniteMPS(4, 45)
+    @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1e-10
+    ψ, envs, δ = find_groundstate(ψ, H, alg)
+    @test E₀(4) ≈ expectation_value(ψ, H, envs) atol = 1e-2
+end
+
+@testset "Z4Irrep" begin
+    H = quantum_potts(Z4Irrep; q=4)
+    ψ = InfiniteMPS(Vect[Z4Irrep](i => 1 for i in 0:3), Vect[Z4Irrep](i => 12 for i in 0:3))
+    @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1e-10
+    ψ, envs, δ = find_groundstate(ψ, H, alg)
+    @test E₀(4) ≈ expectation_value(ψ, H, envs) atol = 1e-2
+end
+
+# Q = 5
+@testset "Q=5 Trivial" begin
+    H = quantum_potts(; q=5)
+    ψ = InfiniteMPS(5, 60)
+    @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1e-10
+    ψ, envs, δ = find_groundstate(ψ, H, alg)
+    @test E₀(5) ≈ expectation_value(ψ, H, envs) atol = 1e-2
+end
+
+@testset "ZNIrrep{5}" begin
+    H = quantum_potts(ZNIrrep{5}; q=5)
+    ψ = InfiniteMPS(Vect[ZNIrrep{5}](i => 1 for i in 0:4),
+                    Vect[ZNIrrep{5}](i => 12 for i in 0:4))
+    @test imag(expectation_value(ψ, H)) ≈ 0 atol = 1e-10
+    ψ, envs, δ = find_groundstate(ψ, H, alg)
+    @test E₀(5) ≈ expectation_value(ψ, H, envs) atol = 1e-2
 end

test/spinoperators.jl

@@ -1,6 +1,6 @@
 using TensorKit
 using TensorOperations
-using LinearAlgebra: tr
+using LinearAlgebra: tr, I
 using TestExtras
 
 ## No symmetry ##
@@ -132,3 +132,47 @@ end
     @test S_exchange(U1Irrep; spin=spin) ≈
           S_xx(U1Irrep; spin=spin) + S_yy(U1Irrep; spin=spin) + S_zz(U1Irrep; spin=spin) rtol = 1e-3
 end
+
+# potts_ZZ test?
+@testset "non-symmetric Q-state potts operators" for Q in 3:5
+    # inferrability
+    X = @inferred potts_X(; q=Q)
+    Z = @inferred potts_Z(; q=Q)
+    ZZ = @inferred potts_ZZ(; q=Q)
+
+    # clock properties
+    @test convert(Array, X^Q) ≈ I
+    @test convert(Array, Z^Q) ≈ I
+
+    # dagger should be reversing the clock direction
+    for s in [X Z]
+        for i in 1:Q
+            @test (s')^i ≈ s^(Q - i)
+        end
+    end
+
+    # commutation relations
+    ω = cis(2π / Q)
+    @test Z * X ≈ ω * X * Z
+end
+
+# potts_ZZ test?
+@testset "Z_Q-symmetric Q-state Potts operators" for Q in 3:5
+    # array conversion
+    _, _, W = weyl_heisenberg_matrices(Q, ComplexF64)
+    @test W * convert(Array, potts_X(; q=Q)) * W' ≈ convert(Array, potts_X(ZNIrrep{Q}; q=Q))
+
+    # inferrability
+    X = @inferred potts_X(ZNIrrep{Q}; q=Q)
+    ZZ = @inferred potts_ZZ(ZNIrrep{Q}; q=Q)
+
+    # unitarity
+    @test X * X' ≈ X' * X
+    @test convert(Array, X * X') ≈ I
+    @test convert(Array, X^Q) ≈ I
+
+    # dagger should be reversing the clock direction
+    for i in 1:Q
+        @test (X')^i ≈ X^(Q - i)
+    end
+end