Benchmark 1

We study the Hamiltonian of the Heisenberg model with periodic boundary conditions.

using CondensedMatterSOS
import MultivariatePolynomials as MP
@spin σ[1:3]
MP.monomials(σ[1], 0:2, consecutive=true)
heisenberg_hamiltonian(σ, true)

\[ σᶻ_{2}σᶻ_{3} + σʸ_{2}σʸ_{3} + σˣ_{2}σˣ_{3} + σᶻ_{1}σᶻ_{3} + σᶻ_{1}σᶻ_{2} + σʸ_{1}σʸ_{3} + σʸ_{1}σʸ_{2} + σˣ_{1}σˣ_{3} + σˣ_{1}σˣ_{2} \]

Let's pick a solver from this list.

import Clarabel
solver = Clarabel.Optimizer

Clarabel.MOIwrapper.Optimizer

We can compute a lower bound -2√2 to the ground state energy as follow:

function hamiltonian_energy(N, maxdegree, solver; symmetry=true, consecutive=false, kws...)
    @spin σ[1:N]
    G = Lattice1Group(N)
    @assert iseven(maxdegree)
    H = heisenberg_hamiltonian(σ, true)
    cone = NonnegPolyInnerCone{MOI.HermitianPositiveSemidefiniteConeTriangle}()
    certificate = SumOfSquares.Certificate.FixedBasis(
        cone,
        MonomialBasis(MP.monomials(σ[1], 0:div(maxdegree, 2), consecutive=consecutive)),
    )
    if symmetry
        certificate = Symmetry.Ideal(
            Symmetry.Pattern(G, Action(σ)),
            certificate,
        )
    end
    energy(H, maxdegree, solver; certificate, kws...)
end
bound, gram, ν = hamiltonian_energy(
    2,
    2,
    solver,
    symmetry = false,
    sparsity = SumOfSquares.Sparsity.NoPattern(),
)
bound

-5.999999993853594

We can see that the moment matrix uses all monomials:

ν.basis.monomials

7-element Vector{CondensedMatterSOS.SpinMonomial}:
 1
 σᶻ₂
 σʸ₂
 σˣ₂
 σᶻ₁
 σʸ₁
 σˣ₁

Symmetry reduction

We can reduce the computation using symmetry reduction as follows.

using CondensedMatterSOS

bound, gram, ν = hamiltonian_energy(
    2,
    2,
    solver,
)
bound

-5.999999997057931

The reduction is obtained by block diagonalizing with a change of polynomial basis to the isotypical basis.

[M.basis.polynomials for M in ν.blocks]

7-element Vector{Vector{MultivariatePolynomials.VectorPolynomial{ComplexF64, MultivariatePolynomials.Term{ComplexF64, CondensedMatterSOS.SpinMonomial}}}}:
 [(-0.7071067811865472 + 0.0im)σᶻ₂ + (-0.7071067811865475 + 0.0im)σᶻ₁]
 [(-0.7071067811865472 + 0.0im)σʸ₂ + (-0.7071067811865475 + 0.0im)σʸ₁]
 [(-0.7071067811865472 + 0.0im)σˣ₂ + (-0.7071067811865475 + 0.0im)σˣ₁]
 [(1.0 - 0.0im)]
 [(-0.7071067811865472 + 0.0im)σᶻ₂ + (0.7071067811865475 + 0.0im)σᶻ₁]
 [(-0.7071067811865472 + 0.0im)σʸ₂ + (0.7071067811865475 + 0.0im)σʸ₁]
 [(-0.7071067811865472 + 0.0im)σˣ₂ + (0.7071067811865475 + 0.0im)σˣ₁]

Let's try this for 3 sites. First without symmetry.

bound, gram, ν = hamiltonian_energy(
    3,
    2,
    solver,
    symmetry = false,
)
@show bound

-4.499999995399682

Now with symmetry.

bound, gram, ν = hamiltonian_energy(
    3,
    2,
    solver,
)
@show bound

-4.499999996767557

Let's look at the isotypical basis.

[M.basis.polynomials for M in ν.blocks]

10-element Vector{Vector{MultivariatePolynomials.VectorPolynomial{ComplexF64, MultivariatePolynomials.Term{ComplexF64, CondensedMatterSOS.SpinMonomial}}}}:
 [(-0.5773502691896257 + 0.0im)σᶻ₃ + (-0.5773502691896257 + 0.0im)σᶻ₂ + (-0.5773502691896257 + 0.0im)σᶻ₁]
 [(-0.5773502691896257 + 0.0im)σʸ₃ + (-0.5773502691896257 + 0.0im)σʸ₂ + (-0.5773502691896257 + 0.0im)σʸ₁]
 [(-0.5773502691896257 + 0.0im)σˣ₃ + (-0.5773502691896257 + 0.0im)σˣ₂ + (-0.5773502691896257 + 0.0im)σˣ₁]
 [(1.0 - 0.0im)]
 [(-0.577350269189626 + 0.0im)σᶻ₃ + (0.2886751345948131 - 0.49999999999999994im)σᶻ₂ + (0.28867513459481275 + 0.5000000000000001im)σᶻ₁]
 [(-0.577350269189626 + 0.0im)σʸ₃ + (0.2886751345948131 - 0.49999999999999994im)σʸ₂ + (0.28867513459481275 + 0.5000000000000001im)σʸ₁]
 [(-0.577350269189626 + 0.0im)σˣ₃ + (0.2886751345948131 - 0.49999999999999994im)σˣ₂ + (0.28867513459481275 + 0.5000000000000001im)σˣ₁]
 [(-0.577350269189626 + 0.0im)σᶻ₃ + (0.28867513459481275 + 0.5000000000000001im)σᶻ₂ + (0.2886751345948131 - 0.49999999999999994im)σᶻ₁]
 [(-0.577350269189626 + 0.0im)σʸ₃ + (0.28867513459481275 + 0.5000000000000001im)σʸ₂ + (0.2886751345948131 - 0.49999999999999994im)σʸ₁]
 [(-0.577350269189626 + 0.0im)σˣ₃ + (0.28867513459481275 + 0.5000000000000001im)σˣ₂ + (0.2886751345948131 - 0.49999999999999994im)σˣ₁]

Now let's define a function for our common use case.

function f(N, d=1; verbose = 1, kws...)
    @show N
    println("***")
    @show d
    bound, _, ν = @time hamiltonian_energy(N, 2d, solver; kws...)
    @show bound
    block_sizes = map(ν.blocks) do M
        return length(M.basis.polynomials)
    end
    @show block_sizes
    if verbose >= 2
        for M in ν.blocks
            display(round.(M.basis.polynomials, digits=6))
        end
    end
    println("E/N = ", bound / N)
    println("------------------------------------")
    return bound / N
end

f (generic function with 2 methods)

With d = 1, we find a lower bound of -3:

lb = f(6, 1, consecutive = true, symmetry = true)

-2.999999998927544

Now with d = 2, we find -2:

lb = f(6, 2, consecutive = true, symmetry = true)

-1.9999999494851501

Now with d = 3, we find -1.8685:

lb = f(6, 3, consecutive = true, symmetry = true)

-1.8495207022514615

This page was generated using Literate.jl.