Anisotropic 2D Ising model

Project.toml

@@ -14,6 +14,7 @@ MatrixAlgebraKit = "6c742aac-3347-4629-af66-fc926824e5e4"
 OptimKit = "77e91f04-9b3b-57a6-a776-40b61faaebe0"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
+Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec"
@@ -29,6 +30,7 @@ LoggingExtras = "~1.0"
 MatrixAlgebraKit = "0.6.1"
 OptimKit = "0.4.0"
 QuadGK = "2.11.2"
+Roots = "3.0.0"
 SpecialFunctions = "2.6.1"
 StableRNGs = "1.0.4"
 TensorKit = "0.16.2"

src/TNRKit.jl

@@ -10,6 +10,7 @@ using DocStringExtensions
 using SpecialFunctions
 using FastGaussQuadrature
 using QuadGK
+using Roots
 using Base.Threads
 using Combinatorics: permutations
 
@@ -93,7 +94,8 @@ include("models/ising.jl")
 include("models/ising_triangular.jl")
 include("models/ising_honeycomb.jl")
 export classical_ising, ising_βc, f_onsager, ising_cft_exact,
-    ising_βc_3D, classical_ising_3D, ising_3D_free_energy_htse, classical_ising_impurity,
+    ising_βc_3D, classical_ising_3D, ising_3D_free_energy_htse,
+    classical_ising_impurity, ising_anisotropic_βc, f_onsager_anisotropic,
     classical_ising_triangular, ising_βc_triangular, f_onsager_triangular,
     classical_ising_honeycomb, ising_βc_honeycomb, f_onsager_honeycomb
 

src/models/ising.jl

@@ -57,76 +57,150 @@ function ising_3D_free_energy_htse(β::Real; J::Real = 1.0, max_order::Int = 24)
 end
 ising_3d_free_energy_htse(; kwargs...) = ising_3D_free_energy_htse(ising_βc_3D; kwargs...)
 
+"""
+    ising_bond_tensor(β::Real, T::Type{<:Number})
+
+Constructs the bond tensor `diag(√cosh(β), √sinh(β))` used to decompose
+the Ising Boltzmann weight on a bond with effective coupling `β`.
+"""
 function ising_bond_tensor(β::Real, T::Type{<:Number})
     x = cosh(β)
     y = sinh(β)
     bond_matrix = T[sqrt(x) 0; 0 sqrt(y)]
     return TensorMap(bond_matrix, ℂ^2 ← ℂ^2)
 end
 
+"""
+    ising_anisotropic_βc(Jx::Real, Jy::Real)
+
+Returns the critical inverse temperature `βc` for the anisotropic 2D Ising model
+on a square lattice, determined by the condition
+
+    sinh(2 βc Jx) · sinh(2 βc Jy) = 1 .
+
+If `Jx == Jy`, this reduces to the isotropic critical point `ising_βc / Jx`.
+"""
+function ising_anisotropic_βc(Jx::Real, Jy::Real)
+    if Jx == Jy
+        return Float64(ising_βc) / Jx
+    end
+    f(β) = sinh(2β * Jx) * sinh(2β * Jy) - 1.0
+    β_max = Float64(ising_βc) / min(Jx, Jy) * 5.0
+    while f(β_max) < 0
+        β_max *= 2.0
+    end
+    return find_zero(f, (0.0, β_max))
+end
+
+"""
+    f_onsager_anisotropic(β::Real, Jx::Real, Jy::Real)
+
+Computes the exact Onsager free energy **per site** for
+the 2D Ising model on a square lattice. The general formula
+with anisotropic couplings `Jx`, `Jy` is
+
+    -βf = ln(2) + 1/(8π²) ∫₀^{2π}∫₀^{2π} dθ₁dθ₂ ln[cosh(2Kx)cosh(2Ky) - sinh(2Kx)cos θ₁ - sinh(2Ky)cos θ₂]
+
+where `Kx = β Jx`, `Ky = β Jy`.
+"""
+function f_onsager_anisotropic(β::Real, Jx::Real, Jy::Real)
+    Kx, Ky = Float64(β * Jx), Float64(β * Jy)
+    # Only use the high-precision constant at the isotropic critical point with J=1
+    if Jx == 1 && Jy == 1 && abs(Kx - Float64(ising_βc)) < 1.0e-14
+        return Float64(f_onsager)
+    end
+
+    c1, s1 = cosh(2Kx), sinh(2Kx)
+    c2, s2 = cosh(2Ky), sinh(2Ky)
+    s2_sq = s2^2
+
+    # The 2D Onsager integral reduces to 1D after integrating out θ₂ analytically
+    # (∫₀^{2π} ln(A - B cos θ) dθ = 2π ln((A + √(A²-B²))/2)):
+    #
+    #   -βf = ln(2) + 1/(2π) ∫₀^{π} dθ  ln((A(θ) + √(A(θ)² - s₂²)) / 2)
+    #
+    # where A(θ) = c₁c₂ - s₁ cos θ.
+    integral, _ = quadgk(0, π) do θ
+        A = c1 * c2 - s1 * cos(θ)
+        log((A + sqrt(A^2 - s2_sq)) / 2)
+    end
+
+    return -(log(2.0) + integral / (2π)) / β
+end
+
 """
     classical_ising(; kwargs...)
     classical_ising(β::Real; kwargs...)
-    classical_ising(::Type{Trivial}, β::Real; T::Type{<:Number} = Float64, h = 0.0)
-    classical_ising(::Type{Z2Irrep}, β::Real; T::Type{<:Number} = Float64, h = 0.0)
+    classical_ising(::Type{Trivial}, β::Real; T::Type{<:Number} = Float64, h = 0.0, Jx = 1.0, Jy = 1.0)
+    classical_ising(::Type{Z2Irrep}, β::Real; T::Type{<:Number} = Float64, h = 0.0, Jx = 1.0, Jy = 1.0)
 
 Constructs the partition function tensor for a 2D square lattice
-for the classical Ising model with a given inverse temperature `β` and external magnetic field `h`.
-Compatible with no symmetry for `h ≠ 0` or with explicit ℤ₂ symmetry for `h = 0` on each of its spaces.
+for the classical Ising model with inverse temperature `β` and external magnetic field `h`.
+Compatible with no symmetry for `h ≠ 0` or with explicit ℤ₂ symmetry for `h = 0`.
 Defaults to ℤ₂ symmetry and `h = 0` if the symmetry type and magnetic field are not provided.
 
-# Examples
+For the anisotropic model, the coupling constants `Jx` (horizontal bonds)
+and `Jy` (vertical bonds) can be specified independently.
+The effective couplings are `Kx = β Jx` and `Ky = β Jy`.
+Defaults to the isotropic case `Jx = Jy = 1.0`.
+
+### Examples
 ```julia
-    classical_ising() # Default symmetry is `Z2Irrep`, default inverse temperature is `ising_βc` and default magnetic field `h = 0`.
-    classical_ising(Trivial, 0.5; h = 1.0) # Custom inverse temperature without symmetry and custom magnetic field `h`.
+    classical_ising()                           # default: ℤ₂ symmetric, isotropic at βc
+    classical_ising(Trivial, 0.5; h = 1.0)      # no symmetry, with magnetic field
+    classical_ising(1.0; Jx = 1.0, Jy = 0.5)    # anisotropic: Jx=1, Jy=0.5
+    classical_ising(Trivial, 0.5; Jy = 0.8)     # anisotropic without symmetry
 
 !!! info
     When studying this model with impurities, the tensor without symmetry should be constructed, as the impurity breaks the ℤ₂ symmetry.
-```
 
-See also: [`classical_ising_3D`](@ref).
+See also: [`classical_ising_3D`](@ref), [`ising_anisotropic_βc`](@ref).
 """
 function classical_ising(β::Real; kwargs...)
     return classical_ising(Z2Irrep, β; kwargs...)
 end
 classical_ising(; kwargs...) = classical_ising(ising_βc; kwargs...)
 classical_ising(::Type{Trivial}; kwargs...) = classical_ising(Trivial, ising_βc; kwargs...)
-function classical_ising(::Type{Trivial}, β::Real; T::Type{<:Number} = Float64, h = 0.0)
+function classical_ising(::Type{Trivial}, β::Real; T::Type{<:Number} = Float64, h = 0.0, Jx = 1.0, Jy = 1.0)
+    Kx, Ky = β * Jx, β * Jy
     init = zeros(T, 2, 2, 2, 2)
     for (i, j, k, l) in Iterators.product([1:2 for _ in 1:4]...)
         init[i, j, k, l] = mod(i + j + k + l, 2) == 0 ? cosh(h * β) : sinh(h * β)
     end
     init = TensorMap(init, ℂ^2 ⊗ ℂ^2 ← ℂ^2 ⊗ ℂ^2)
 
-    bond_tensor = ising_bond_tensor(β, T)
-
-    @tensor T[-1 -2; -3 -4] := 2 * init[1 2; 3 4] * bond_tensor[-1; 1] * bond_tensor[-2; 2] * bond_tensor[3; -3] * bond_tensor[4; -4]
+    bond_tensor_x = ising_bond_tensor(Kx, T)
+    bond_tensor_y = ising_bond_tensor(Ky, T)
+    @tensor T[-1 -2; -3 -4] := 2 * init[1 2; 3 4] * bond_tensor_x[-1; 1] * bond_tensor_y[-2; 2] * bond_tensor_y[3; -3] * bond_tensor_x[4; -4]
     return T
 end
-function classical_ising(::Type{Z2Irrep}, β::Real; T::Type{<:Number} = Float64, h = 0.0)
+function classical_ising(::Type{Z2Irrep}, β::Real; T::Type{<:Number} = Float64, h = 0.0, Jx = 1.0, Jy = 1.0)
     @assert h == 0.0 "External magnetic field is not compatible with ℤ₂ symmetry"
-    x = cosh(β)
-    y = sinh(β)
+    Kx, Ky = β * Jx, β * Jy
+    xh, yh = cosh(Kx), sinh(Kx)
+    xv, yv = cosh(Ky), sinh(Ky)
+    w = sqrt(xh * yh * xv * yv)
 
     S = ℤ₂Space(0 => 1, 1 => 1)
     t = zeros(T, S ⊗ S ← S ⊗ S)
-    block(t, Irrep[ℤ₂](0)) .= [2x^2 2x * y; 2x * y 2y^2]
-    block(t, Irrep[ℤ₂](1)) .= [2x * y 2x * y; 2x * y 2x * y]
+    block(t, Irrep[ℤ₂](0)) .= [2 * xh * xv 2 * w; 2 * w 2 * yh * yv]
+    block(t, Irrep[ℤ₂](1)) .= [2 * w 2 * yh * xv; 2 * xh * yv 2 * w]
 
     return t
 end
 
 """
-    classical_ising_impurity([Type{Trivial}], β::Real; T::Type{<:Number} = Float64, h = 0.0)
+    classical_ising_impurity([Type{Trivial}], β::Real; T::Type{<:Number} = Float64, h = 0.0, Jx = 1.0, Jy = 1.0)
 
 Constructs the partition function tensor for a 2D square lattice
-for the classical Ising model with a given inverse temperature `β` and external magnetic field `h` with a magnetisation impurity.
-Compatible with no symmetry on each of its spaces.
+for the classical Ising model with a given inverse temperature `β` and external magnetic field `h`
+with a magnetisation impurity. Compatible with no symmetry on each of its spaces.
 
 # Examples
 ```julia
-    classical_ising_impurity() # Default inverse temperature is `ising_βc`
-    classical_ising_impurity(0.5; h = 1.0) # Custom inverse temperature and magnetic field
+    classical_ising_impurity()                          # default: isotropic at βc
+    classical_ising_impurity(0.5; h = 1.0)              # with magnetic field
+    classical_ising_impurity(0.5; Jx = 1.0, Jy = 0.5)   # anisotropic couplings
 ```
 !!! info
     When calculating the free energy with `free_energy()`, set the `initial_size` keyword argument to `2.0`.
@@ -138,16 +212,18 @@ function classical_ising_impurity(β::Real; kwargs...)
     return classical_ising_impurity(Trivial, β; kwargs...)
 end
 classical_ising_impurity(; kwargs...) = classical_ising_impurity(ising_βc; kwargs...)
-function classical_ising_impurity(::Type{Trivial}, β::Real; T::Type{<:Number} = Float64, h = 0.0)
+function classical_ising_impurity(::Type{Trivial}, β::Real; T::Type{<:Number} = Float64, h = 0.0, Jx = 1.0, Jy = 1.0)
+    Kx, Ky = β * Jx, β * Jy
     init = zeros(T, 2, 2, 2, 2)
     for (i, j, k, l) in Iterators.product([1:2 for _ in 1:4]...)
         init[i, j, k, l] = mod(i + j + k + l, 2) == 0 ? sinh(h * β) : cosh(h * β)
     end
     init = TensorMap(init, ℂ^2 ⊗ ℂ^2 ← ℂ^2 ⊗ ℂ^2)
 
-    bond_tensor = ising_bond_tensor(β, T)
+    bond_tensor_x = ising_bond_tensor(Kx, T)
+    bond_tensor_y = ising_bond_tensor(Ky, T)
 
-    @tensor t[-1 -2; -3 -4] := 2 * init[1 2; 3 4] * bond_tensor[-1; 1] * bond_tensor[-2; 2] * bond_tensor[3; -3] * bond_tensor[4; -4]
+    @tensor t[-1 -2; -3 -4] := 2 * init[1 2; 3 4] * bond_tensor_x[-1; 1] * bond_tensor_y[-2; 2] * bond_tensor_y[3; -3] * bond_tensor_x[4; -4]
     return t
 end
 

test/models/models.jl

@@ -7,9 +7,15 @@ println("--------------------")
 println(" Testing all models ")
 println("--------------------")
 
+# Pre-computed anisotropic Ising critical values for testing (Jx=1, Jy=0.5)
+const ising_aniso_βc_1_05 = 0.6093778634360062
+const f_onsager_aniso_1_05 = -1.5741477440817498
+
 model_temp_answer_string_2d = [
     (classical_ising(Trivial), ising_βc, f_onsager, "2D Ising model with no symmetry"),
     (classical_ising(), ising_βc, f_onsager, "2D Ising model with ℤ₂ symmetry"),
+    (classical_ising(ising_aniso_βc_1_05; Jx = 1.0, Jy = 0.5), ising_aniso_βc_1_05, f_onsager_aniso_1_05, "2D anisotropic Ising model Jx=1,Jy=0.5 with ℤ₂ symmetry"),
+    (classical_ising(Trivial, ising_aniso_βc_1_05; Jx = 1.0, Jy = 0.5), ising_aniso_βc_1_05, f_onsager_aniso_1_05, "2D anisotropic Ising model Jx=1,Jy=0.5 with no symmetry"),
     (gross_neveu_start(0, 0, 0), 1.0, -1.4515448845652446, "Gross-Neveu model"),
     (classical_clock(Trivial, 3, 2.0 * log(√3 + 1) / 3), 2.0 * log(√3 + 1) / 3, -4.17924244901635, "3-state clock model with no symmetry"), # This is an approximation!
     (classical_clock(Z3Irrep, 3, 2.0 * log(√3 + 1) / 3), 2.0 * log(√3 + 1) / 3, -4.17924244901635, "3-state clock model with ℤ₃ symmetry"), # This is an approximation!
@@ -82,6 +88,33 @@ end
     @info "Obtained central charge:\n$central_charge."
 end
 
+# Tests for anisotropic Ising helper functions
+@testset "Anisotropic Ising helpers" begin
+    # Critical condition: sinh(2βc Jx) · sinh(2βc Jy) = 1, and βc scaling
+    for (Jx, Jy) in [(1.0, 1.0), (1.0, 0.5), (1.0, 0.3), (2.0, 2.0)]
+        βc = ising_anisotropic_βc(Jx, Jy)
+        @test sinh(2 * Float64(βc) * Jx) * sinh(2 * Float64(βc) * Jy) ≈ 1.0 rtol = 1.0e-14
+        if Jx == Jy
+            @test Float64(βc) ≈ Float64(ising_βc) / Jx  # βc scales as 1/J
+        end
+    end
+
+    # Free energy scaling: f(β, J, J) = J · f_iso(β·J)
+    for (β, J) in [(Float64(ising_βc), 1.0), (0.5, 2.0), (0.3, 1.5)]
+        βc = ising_anisotropic_βc(J, J)
+        f_J = f_onsager_anisotropic(β, J, J)
+        f_Jc = f_onsager_anisotropic(βc, J, J)
+        f_1 = f_onsager_anisotropic(β * J, 1.0, 1.0)
+        @test f_J ≈ J * f_1 rtol = 1.0e-6
+        @test f_Jc ≈ J * Float64(f_onsager) rtol = 1.0e-6
+    end
+
+    # Jy→0 limit: f should approach the 1D Ising result
+    β = 0.5
+    f_1d = -(log(2.0) + log(cosh(β))) / β
+    @test f_onsager_anisotropic(β, 1.0, 1.0e-10) ≈ f_1d rtol = 1.0e-6
+end
+
 for (model, temp, answer, description) in model_temp_answer_string_3d
     @testset "$(description)" begin
         scheme = HOTRG_3D(model)

test/runtests.jl

@@ -1,6 +1,12 @@
-using TNRKit
 using ParallelTestRunner
+using TNRKit
+
+# --fast to indicate a smaller set of tests
+args = parse_args(ARGS; custom = ["fast"])
+fast = !isnothing(args.custom["fast"])
+
+const init_code = quote
+    const fast_tests = $fast
+end
 
-testsuite = find_tests(@__DIR__)
-args = parse_args(ARGS)
-ParallelTestRunner.runtests(TNRKit, args)
+ParallelTestRunner.runtests(TNRKit, args; init_code)