From f474691926e762eccbd013cee3b91fca569cf85e Mon Sep 17 00:00:00 2001 From: Adwait-Naravane Date: Thu, 10 Sep 2026 16:50:14 +0200 Subject: [PATCH] Add Initial tensor for Complex scalar field theory in 2D with O(2) symmetry --- README.md | 2 +- docs/src/index.md | 2 +- src/models/phi4_complex.jl | 167 ++++++++++++++++++++++++++++++++----- test/models/models.jl | 32 +++++++ 4 files changed, 182 insertions(+), 21 deletions(-) diff --git a/README.md b/README.md index 47891018..8b2cd735 100644 --- a/README.md +++ b/README.md @@ -119,7 +119,7 @@ TNRKit includes several common models out of the box. - XY model in 2D: `classical_XY(S, β, charge_trunc)` where `S` can be `U1Irrep` or `CU1Irrep` to specify the symmetry. - Real $\phi^4$ model: `phi4_real(S, K, μ0, λ, h)` where `S` can be `Trivial` or `Z2Irrep` to specify the symmetry. - Real $\phi^4$ model with impurities: `phi4_real_imp1(S, K, μ0, λ, h)` and `phi4_real_imp2(S, K, μ0, λ, h)` where `S` can be `Trivial`. -- Complex $\phi^4$ model: `phi4_complex(S, K, μ0, λ)` where `S` can be `Trivial`, `Z2Irrep ⊠ Z2Irrep` or `U1Irrep` to specify the symmetry. +- Complex $\phi^4$ model: `phi4_complex(S, K, μ0, λ)` where `S` can be `Trivial`, `Z2Irrep ⊠ Z2Irrep`, `U1Irrep` or `CU1Irrep` to specify the symmetry. - Gross-Neveu model: `gross_neveu_start(S, μ, m, g)` where `S` can be `FermionParity` to specify the symmetry. ## Included Models on the triangular lattice diff --git a/docs/src/index.md b/docs/src/index.md index 0cd019e5..24ef0ee6 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -93,7 +93,7 @@ TNRKit includes several common models out of the box. - XY model in 2D: `classical_XY(S, β, charge_trunc)` where `S` can be `U1Irrep` or `CU1Irrep` to specify the symmetry. - Real $\phi^4$ model: `phi4_real(S, K, μ0, λ, h)` where `S` can be `Trivial` or `Z2Irrep` to specify the symmetry. - Real $\phi^4$ model with impurities: `phi4_real_imp1(S, K, μ0, λ, h)` and `phi4_real_imp2(S, K, μ0, λ, h)` where `S` can be `Trivial`. -- Complex $\phi^4$ model: `phi4_complex(S, K, μ0, λ)` where `S` can be `Trivial`, `Z2Irrep ⊠ Z2Irrep` or `U1Irrep` to specify the symmetry. +- Complex $\phi^4$ model: `phi4_complex(S, K, μ0, λ)` where `S` can be `Trivial`, `Z2Irrep ⊠ Z2Irrep`, `U1Irrep` or `CU1Irrep` to specify the symmetry. - Gross-Neveu model: `gross_neveu_start(S, μ, m, g)` where `S` can be `FermionParity` to specify the symmetry. ## Included Models on the triangular lattice diff --git a/src/models/phi4_complex.jl b/src/models/phi4_complex.jl index 6ac7937a..2b1def9c 100644 --- a/src/models/phi4_complex.jl +++ b/src/models/phi4_complex.jl @@ -6,8 +6,10 @@ function f_complex(ℝϕ1::Float64, ℂϕ1::Float64, ℝϕ2::Float64, ℂϕ2::Float64, μ0::Float64, λ::Float64) return exp( -1 / 2 * ((ℝϕ1 - ℝϕ2)^2 + (ℂϕ1 - ℂϕ2)^2) - - μ0 / 8 * (ℝϕ1^2 + ℂϕ1^2 + ℝϕ2^2 + ℂϕ2^2) - - λ / 16 * ((ℝϕ1^2 + ℂϕ1^2)^2 + (ℝϕ2^2 + ℂϕ2^2)^2) + - + μ0 / 8 * (ℝϕ1^2 + ℂϕ1^2 + ℝϕ2^2 + ℂϕ2^2) + - + λ / 16 * ((ℝϕ1^2 + ℂϕ1^2)^2 + (ℝϕ2^2 + ℂϕ2^2)^2) ) end @@ -54,6 +56,60 @@ function precompute_moments_complex(K, μ0, λ) return M end +# For phi4_complex_U1 and phi4_complex_CU1 +# `logfact[n + 1] == log(n!)`, accumulated in log space so that it stays finite +# for the large `K` where `factorial(n)` would overflow. +function phi4_complex_logfactorials(K) + return [0.0; cumsum(log.(1.0:(K - 1)))] +end + +# For phi4_complex_U1 and phi4_complex_CU1 +# A single entry of the Taylor-expanded tensor in the exponent basis, in which a +# leg state |a, b⟩ carries `a` powers of ϕ and `b` powers of ϕ̄: +# +# 2π δ(a + c + f + h - b - d - e - g) +# ─────────────────────────────────────────────── × ∫₀^∞ dr r^{a+b+…+h+1} +# √(2^{a+b+…+h} a! b! c! d! e! f! g! h!) × exp(-(2 + µ_0²/2) r² - (λ/4) r⁴) +# +# with the radial integral supplied by `precompute_moments_complex`. +function phi4_complex_weight(moments, logfact, a, b, c, d, e, f, g, h) + # The δ, i.e. U(1) charge conservation: (a - b) + (c - d) == (e - f) + (g - h) + (a + c + f + h == b + d + e + g) || return 0.0 + + n = a + b + c + d + e + f + g + h + M = moments[n + 2] + iszero(M) && return 0.0 + + logdenom = 0.5 * ( + log(2) * n + + logfact[a + 1] + logfact[b + 1] + logfact[c + 1] + logfact[d + 1] + + logfact[e + 1] + logfact[f + 1] + logfact[g + 1] + logfact[h + 1] + ) + return 2π * M / exp(logdenom) +end + +# For phi4_complex_CU1 +# CU(1) = O(2) adapted leg space. Charge conjugation acts on the exponent basis +# as C|a, b⟩ = |b, a⟩, so the states organise into irreps as +# * a - b == 0: |a, a⟩ is C-even -> sector (0, 0), multiplicity K +# * a - b == q > 0: {|b+q, b⟩, |b, b+q⟩} -> sector (q, 2), multiplicity K - q +# The C-odd sector (0, 1) does not occur. The total dimension is +# K + 2 Σ_{q=1}^{K-1} (K - q) = K^2, the same bond dimension as the U(1) tensor. +function phi4_complex_cu1_space(K) + return CU1Space(vcat([(0, 0) => K], [(q, 2) => K - q for q in 1:(K - 1)])) +end + +# For phi4_complex_CU1 +# Exponents (a, b) of the leg state at position `i` inside CU1 sector `s` with +# multiplicity index `m`. TensorKitSectors fixes the basis of a two-dimensional +# (j, 2) irrep to be (index 1, index 2) = (charge +j, charge -j) -- see the +# `fusiontensor(::CU1Irrep, ...)` branch `c.j == a.j + b.j` +function phi4_complex_cu1_exponents(s::CU1Irrep, i::Integer, m::Integer) + q = convert(Int, s.j) + q == 0 && return (m - 1, m - 1) + return i == 1 ? (m - 1 + q, m - 1) : (m - 1, m - 1 + q) +end + # For phi4_complex_Z2Z2 function precompute_radial_integrals(N, μ0, λ; rtol = 1.0e-8) @@ -147,12 +203,13 @@ end phi4_complex(::Type{Trivial}, K::Integer, μ0::Float64, λ::Float64; T::Type{<:Number} = Float64) phi4_complex(::Type{Z2Irrep ⊠ Z2Irrep}, K::Integer, μ0::Float64, λ::Float64; T::Type{<:Number} = Float64) phi4_complex(::Type{U1Irrep}, K::Integer, μ0::Float64, λ::Float64; T::Type{<:Number} = Float64) + phi4_complex(::Type{CU1Irrep}, K::Integer, μ0::Float64, λ::Float64; T::Type{<:Number} = Float64) Constructs the partition function tensor for a 2D square lattice for the complex ϕ^4 model with a given approximation `K`, bare mass µ_0^2 `μ0` and interaction constant `λ`. It is based on [Gauss-Hermite quadrature](https://en.wikipedia.org/wiki/Gauss%E2%80%93Hermite_quadrature). -Compatible with no symmetry, explicit ℤ₂×ℤ₂ symmetry or explicit U(1) symmetry on each of its spaces. +Compatible with no symmetry, explicit ℤ₂×ℤ₂ symmetry, explicit U(1) symmetry or explicit CU(1) = O(2) symmetry on each of its spaces. Defaults to U(1) symmetry if the symmetry type is not provided. # Arguments @@ -173,6 +230,21 @@ The order of the Taylor expansion is `K`. The total bond dimension is `K^2`. The tensor is constructed by Taylor expanding the mixed sites term in the partition function. The order of the Taylor expansion is `K`. The total bond dimension is `K^2`. +Every leg carries a pair of Taylor exponents `(a, b)`, it counts the powers of ϕ and of ϕ̄, so +it carries the U(1) charge `q = a - b`. + +## CU(1) = O(2) symmetry +Charge conjugation `C: ϕ ↔ ϕ̄` is a symmetry of the model which combined with U(1) gives the O(2) symmetry group. +The tensor is invariant under the charge-conjugation operation, which is just the global swap + `(a, c, e, g) ↔ (b, d, f, h)`, so the U(1) tensor is in fact O(2) = U(1) ⋊ C +symmetric. + +Since `C|a, b⟩ = |b, a⟩` on a leg, the leg space is +`V = (0, 0)^K ⊕ (1, 2)^{K-1} ⊕ … ⊕ (K-1, 2)^1`: the neutral states `|a, a⟩` are C-even and +the charge `±q` states pair up into the two-dimensional irrep `(q, 2)`. The C-odd sector +`(0, 1)` does not appear. The total bond dimension is again `K^2`, but each `±q` pair of +U(1) blocks is merged into a single block, roughly halving the amount of stored data. + # Examples ```julia phi4_complex(10, -1., 1.) @@ -182,7 +254,7 @@ The order of the Taylor expansion is `K`. The total bond dimension is `K^2`. When studying this model with impurities, the tensor without symmetry should be constructed, as the impurity breaks the symmetry. # References -Piceu Jarid and Adwait Naravane, but based on: +Jarid Piceu and Adwait Naravane, but based on: * [Kadoh et. al. 10.1007/JHEP05(2019)184 (2019)](@cite kadoh2019) * [Delcamp et. al. Phys. Rev. Research 2, 033278 (2020)](@cite delcamp2020) @@ -270,7 +342,7 @@ function phi4_complex(::Type{U1Irrep}, K::Integer, μ0::Float64, λ::Float64; T: end moments = precompute_moments_complex(K, μ0, λ) - logfact = log.(factorial.(0:(K - 1))) + logfact = phi4_complex_logfactorials(K) V1 = U1Space([U1Irrep(q) => 1 for q in 0:(K - 1)]...) V2 = U1Space([U1Irrep(q) => 1 for q in 0:-1:(-K + 1)]...) @@ -303,32 +375,28 @@ function phi4_complex(::Type{U1Irrep}, K::Integer, μ0::Float64, λ::Float64; T: # index as block[i1, i2, i3, i4] where each i is the multiplicity index for a in max(0, q1):min(K - 1, q1 + K - 1) - B = a - q1; (0 <= B <= K - 1) || continue + B = a - q1 + (0 <= B <= K - 1) || continue i1 = mult_index[q1][a] for c in max(0, q2):min(K - 1, q2 + K - 1) - D = c - q2; (0 <= D <= K - 1) || continue + D = c - q2 + (0 <= D <= K - 1) || continue i2 = mult_index[q2][c] for e in max(0, q3):min(K - 1, q3 + K - 1) - F = e - q3; (0 <= F <= K - 1) || continue + F = e - q3 + (0 <= F <= K - 1) || continue i3 = mult_index[q3][e] for g in max(0, q4):min(K - 1, q4 + K - 1) - H = g - q4; (0 <= H <= K - 1) || continue + H = g - q4 + (0 <= H <= K - 1) || continue i4 = mult_index[q4][g] - sum_power = a + B + c + D + e + F + g + H - M = moments[sum_power + 2] - (M == 0.0) && continue - - logdenom = 0.5 * ( - log(2) * sum_power + - logfact[a + 1] + logfact[B + 1] + logfact[c + 1] + logfact[D + 1] + - logfact[e + 1] + logfact[F + 1] + logfact[g + 1] + logfact[H + 1] + block[i1, i2, i3, i4] += phi4_complex_weight( + moments, logfact, a, B, c, D, e, F, g, H ) - - block[i1, i2, i3, i4] += 2π * M / exp(logdenom) end end end @@ -338,6 +406,67 @@ function phi4_complex(::Type{U1Irrep}, K::Integer, μ0::Float64, λ::Float64; T: return T_fused end +function phi4_complex(::Type{CU1Irrep}, K::Integer, μ0::Float64, λ::Float64; T::Type{<:Number} = Float64) + if K % 2 != 0 + error("K must be even") + end + + moments = precompute_moments_complex(K, μ0, λ) + logfact = phi4_complex_logfactorials(K) + + V = phi4_complex_cu1_space(K) + t = zeros(T, V ⊗ V ← V ⊗ V) + + trees = collect(fusiontrees(t)) + @threads for (split_tree, fuse_tree) in trees + s1, s2 = split_tree.uncoupled + s3, s4 = fuse_tree.uncoupled + c = split_tree.coupled + + # Clebsch-Gordan tensors of the two vertices. None of the legs is dual, + # so these are the bare `fusiontensor`s of the sectors involved. + C12 = fusiontensor(s1, s2, c) + C34 = fusiontensor(s3, s4, c) + + block = t[split_tree, fuse_tree] + # block has shape (mult(s1) × mult(s2)) × (mult(s3) × mult(s4)) + + for i1 in 1:dim(s1), i2 in 1:dim(s2), i3 in 1:dim(s3), i4 in 1:dim(s4) + # The projector onto this pair of fusion trees, normalised as in + # `TensorKit.project_symmetric!`. + w = zero(eltype(C12)) + for μ in 1:dim(c) + w += C12[i1, i2, μ, 1] * C34[i3, i4, μ, 1] + end + iszero(w) && continue + w /= dim(c) + + for m1 in axes(block, 1) + a, b = phi4_complex_cu1_exponents(s1, i1, m1) + + for m2 in axes(block, 2) + cc, d = phi4_complex_cu1_exponents(s2, i2, m2) + + for m3 in axes(block, 3) + e, f = phi4_complex_cu1_exponents(s3, i3, m3) + + for m4 in axes(block, 4) + g, h = phi4_complex_cu1_exponents(s4, i4, m4) + + v = phi4_complex_weight(moments, logfact, a, b, cc, d, e, f, g, h) + iszero(v) && continue + + block[m1, m2, m3, m4] += w * v + end + end + end + end + end + end + + return t +end + """ phi4_complex_impϕ([Type{Trivial}], K::Integer, μ0::Float64, λ::Float64; T::Type{<:Number} = ComplexF64) diff --git a/test/models/models.jl b/test/models/models.jl index d4ccbbd5..d089483d 100644 --- a/test/models/models.jl +++ b/test/models/models.jl @@ -2,6 +2,7 @@ using Test using TNRKit using TensorKit using TensorKitSectors +using LinearAlgebra println("--------------------") println(" Testing all models ") @@ -35,6 +36,8 @@ model_temp_answer_string_2d = [ (phi4_complex(Trivial, 6, -1.0, 1.0), -1.0, 0.7583605364656325, "Complex φ⁴ model with no symmetry"), # This is an approximation! (phi4_complex(6, -1.0, 1.0), -1.0, 0.7673189874157453, "Complex φ⁴ model with U(1) symmetry"), # This is an approximation! (phi4_complex(Z2Irrep ⊠ Z2Irrep, 6, -1.0, 1.0), -1.0, 0.7665677554973079, "Complex φ⁴ model with ℤ₂ × ℤ₂ symmetry"), # This is an approximation! + (phi4_complex(U1Irrep, 6, -1.0, 1.0), -1.0, 0.7673189874157453, "Complex φ⁴ model with U(1) symmetry"), # This is an approximation! + (phi4_complex(CU1Irrep, 6, -1.0, 1.0), -1.0, 0.7673190424140406, "Complex φ⁴ model with CU(1) symmetry"), # This is an approximation! ] model_temp_answer_string_3d = [ @@ -50,6 +53,35 @@ for (model, temp, answer, description) in model_temp_answer_string_2d end end +@testset "Complex φ⁴ - CU(1) is the U(1) tensor in an O(2) adapted basis" begin + K, μ0, λ = 4, -1.0, 1.0 + T_u1 = phi4_complex(U1Irrep, K, μ0, λ) + T_cu1 = phi4_complex(CU1Irrep, K, μ0, λ) + + # Same bond dimension, but every ±q pair of U(1) blocks is merged into one. + @test dim(space(T_cu1, 1)) == K^2 + @test dim(space(T_cu1, 1)) == dim(space(T_u1, 1)) + @test length(blocks(T_cu1)) < length(blocks(T_u1)) + + # Projecting onto CU(1) discards whatever is not O(2) symmetric, so an + # unchanged norm is exactly the statement that nothing was thrown away. + @test norm(T_cu1) ≈ norm(T_u1) + + # Basis independent contractions: the 1×1 and 2×2 tori. + @tensor z_u1 = T_u1[a b; a b] + @tensor z_cu1 = T_cu1[a b; a b] + @test z_cu1 ≈ z_u1 + + @tensor z2_u1 = T_u1[a b; c d] * T_u1[c e; a f] * T_u1[g f; h b] * T_u1[h d; g e] + @tensor z2_cu1 = T_cu1[a b; c d] * T_cu1[c e; a f] * T_cu1[g f; h b] * T_cu1[h d; g e] + @test z2_cu1 ≈ z2_u1 + + # ... and the full spectrum of T seen as a (V ⊗ V) × (V ⊗ V) matrix. + σ_u1 = sort(svdvals(reshape(convert(Array, T_u1), K^4, K^4)); rev = true) + σ_cu1 = sort(svdvals(reshape(convert(Array, T_cu1), K^4, K^4)); rev = true) + @test σ_cu1 ≈ σ_u1 +end + @testset "LoopTNR - 2D XY model" begin @info "Central charge of KT phase with U(1) symmetry" T_KT = classical_XY(U1Irrep, XY_βc + 0.1, 8)