diff --git a/docs/src/examples/excitations/su3-heisenberg/figure-1.png b/docs/src/examples/excitations/su3-heisenberg/figure-1.png new file mode 100644 index 000000000..0b7b9b930 Binary files /dev/null and b/docs/src/examples/excitations/su3-heisenberg/figure-1.png differ diff --git a/docs/src/examples/excitations/su3-heisenberg/figure-2.png b/docs/src/examples/excitations/su3-heisenberg/figure-2.png new file mode 100644 index 000000000..ef40ea9ff Binary files /dev/null and b/docs/src/examples/excitations/su3-heisenberg/figure-2.png differ diff --git a/docs/src/examples/excitations/su3-heisenberg/index.md b/docs/src/examples/excitations/su3-heisenberg/index.md new file mode 100644 index 000000000..a9f10a0b8 --- /dev/null +++ b/docs/src/examples/excitations/su3-heisenberg/index.md @@ -0,0 +1,237 @@ +```@meta +EditURL = "../../../../../examples/excitations/su3-heisenberg/main.jl" +``` + +[![](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/QuantumKitHub/MPSKit.jl/gh-pages?filepath=dev/examples/excitations/su3-heisenberg/main.ipynb) +[![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](https://nbviewer.jupyter.org/github/QuantumKitHub/MPSKit.jl/blob/gh-pages/dev/examples/excitations/su3-heisenberg/main.ipynb) +[![](https://img.shields.io/badge/download-project-orange)](https://minhaskamal.github.io/DownGit/#/home?url=https://github.com/QuantumKitHub/MPSKit.jl/examples/tree/gh-pages/dev/examples/excitations/su3-heisenberg) + +# The SU(3) Heisenberg chain + +The spin-1 Heisenberg chain has a gapped, symmetry-protected topological ground state, the Haldane phase. +In this example we study what happens when the SU(2) symmetry of the spins is promoted to SU(3), following [devos2022](@cite). +Every site of the chain carries the fully symmetric ``[3 0 0]`` irrep of SU(3), which is ten-dimensional, and neighbouring sites are coupled by the SU(3)-invariant Heisenberg interaction + +```math +H = J ∑_{⟨i,j⟩} ∑_{a=1}^{8} Tᵃᵢ Tᵃⱼ +``` + +where the ``Tᵃ`` are the eight generators of SU(3). +This is the direct analogue of the spin-1 chain, with the three generators of SU(2) replaced by the eight generators of SU(3). + +Just as in the [Haldane gap](@ref "The Haldane gap") example, imposing the full non-abelian symmetry block-diagonalizes every tensor, so that a modest number of variational parameters captures a large bond dimension. +The SU(3) Clebsch–Gordan coefficients are provided by [SUNRepresentations.jl](https://github.com/QuantumKitHub/SUNRepresentations.jl), whose `SU3Irrep` sectors plug directly into TensorKit. + +!!! warning "First-run cost" + SUNRepresentations.jl computes Clebsch–Gordan coefficients on demand, and caches them both in memory and on disk. + The first time a coefficient is needed, this computation can be slow. + +To follow this example you need the following packages: + +````julia +using MPSKit, TensorKit, SUNRepresentations, CairoMakie +```` + +For reproducibility of this page, we fix the seed of the random number generator: + +````julia +using Random +Random.seed!(4321); +```` + +## The model + +Each site carries the ``[3 0 0]`` irrep, which we specify through its highest weight. + +````julia +I300 = SU3Irrep(3, 0, 0) +Vphys = Vect[SU3Irrep](I300 => 1) +dim(I300) +```` + +```` +10 +```` + +To build the interaction, we use that it is diagonal in the total charge ``c`` of a bond. +Writing ``C₂ = ∑ₐ TᵃTᵃ`` for the quadratic Casimir, we have + +```math +∑ₐ Tᵃᵢ Tᵃⱼ = ½ (C₂(c) - 2 C₂([3 0 0])). +``` + +SUNRepresentations.jl provides the quadratic Casimir as `casimir(2, irrep)`, normalized such that the fundamental ``[1 0 0]`` has ``C₂ = 4/3``. + +````julia +casimir(2, SU3Irrep(1, 0, 0)), casimir(2, I300) +```` + +```` +(4//3, 6//1) +```` + +The nearest-neighbour term is then a single operator that acts as a scalar on each fusion channel of ``[3 0 0] ⊗ [3 0 0]``. +We construct it by filling in the blocks, one for every total charge ``c``: + +````julia +function su3_heisenberg(Vphys; J = 1.0) + c2_site = casimir(2, SU3Irrep(3, 0, 0)) + SS = zeros(ComplexF64, Vphys ⊗ Vphys ← Vphys ⊗ Vphys) + for (c, data) in blocks(SS) + coupling = (casimir(2, c) - 2 * c2_site) / 2 + for i in axes(data, 1) + data[i, i] = J * coupling + end + end + return InfiniteMPOHamiltonian(SS) +end + +H = su3_heisenberg(Vphys) +```` + +```` +1-site InfiniteMPOHamiltonian(ComplexF64, Rep[SU₃]) with maximal dimension 10: +| ⋮ +| (((0, 0, 0)=>1) ⊞ ((2, 1, 0)=>1) ⊞ ((0, 0, 0)=>1)) +┼─[1]─ ((3, 0, 0) => 1) +│ (((0, 0, 0)=>1) ⊞ ((2, 1, 0)=>1) ⊞ ((0, 0, 0)=>1)) +| ⋮ + +```` + +The product ``[3 0 0] ⊗ [3 0 0]`` decomposes into four irreps, and the antiferromagnetic coupling favours the channels with the lowest Casimir: + +````julia +for (c, n) in directproduct(I300, I300) + coupling = (casimir(2, c) - 2 * casimir(2, I300)) / 2 + println("$c:\t dim = $(dim(c)),\t coupling = $coupling") +end +```` + +```` +Irrep[SU₃]((4, 2, 0)): dim = 27, coupling = -2//1 +Irrep[SU₃]((5, 1, 0)): dim = 35, coupling = 0//1 +Irrep[SU₃]((6, 0, 0)): dim = 28, coupling = 3//1 +Irrep[SU₃]((3, 3, 0)): dim = 10, coupling = -3//1 + +```` + +## Ground state + +The ``[3 0 0]`` irrep has trivial triality, so that, as for the spin-1 chain, a translation-invariant state with a single-site unit cell is allowed. +The physical charge then cannot change the triality of the virtual spaces, so all virtual charges must have the same triality. +Here we use virtual charges with trivial triality, which give the lowest variational energy. +We start from a small virtual space, and optimize a random initial state with VUMPS: + +````julia +Vvirt = Vect[SU3Irrep]( + SU3Irrep(0, 0, 0) => 3, + SU3Irrep(2, 1, 0) => 6, + SU3Irrep(3, 0, 0) => 3, + SU3Irrep(3, 3, 0) => 3, + SU3Irrep(4, 2, 0) => 2, +) +ψ₀ = InfiniteMPS([Vphys], [Vvirt]) +vumps = VUMPS(; maxiter = 400, tol = 1.0e-10, verbosity = 0) +ψ, envs, info = find_groundstate(ψ₀, H, vumps) +E₀ = real(expectation_value(ψ, H, envs)) +```` + +```` +-2.5735056524313213 +```` + +## Excitations + +On top of the ground state, we compute excitations with the quasiparticle ansatz. +Following [devos2022](@cite), we look for the lowest excitations in the adjoint ``[2 1 0]`` sector, which has its minimum at momentum ``k = 2π/3``. + +````julia +sector = SU3Irrep(2, 1, 0) +Es, _ = excitations(H, QuasiparticleAnsatz(), 2π / 3, ψ, envs; sector) +ΔE = real(Es[1]) +```` + +```` +-0.04185666026343908 +```` + +The excitation energy comes out negative, which would mean that there are states below the ground state. + +!!! warning "Excitation energies are not variational" + The quasiparticle ansatz returns the energy of the excited state *relative to the ground-state energy*. + The ground state itself is only approximated, and its energy ``E₀`` lies above the exact one. + Any error in ``E₀`` therefore directly shifts the excitation energies, and they are not an upper bound on the true gap. + In particular, a gap that is small compared to the error in ``E₀`` can come out too small, or even negative. + +The remedy is to improve the ground state. +Instead of guessing a larger virtual space, we let `changebonds` with `OptimalExpand` grow the bond dimension, by adding the directions that contribute the most to a two-site update of the converged state. +Here `truncrank` counts the total dimension that is added to each bond, including the dimension of the SU(3) irreps. +Every expanded state is then reconverged with VUMPS, and we recompute the excitation energy at ``k = 2π/3``. + +````julia +expand = OptimalExpand(; trunc = truncrank(150)) +Ds = [dim(left_virtualspace(ψ, 1))] +E₀s = [E₀] +ΔEs = [ΔE] +for step in 1:8 + global ψ, envs = changebonds(ψ, H, expand, envs) + global ψ, envs, info = find_groundstate(ψ, H, vumps, envs) + Es_step, _ = excitations(H, QuasiparticleAnsatz(), 2π / 3, ψ, envs; sector) + push!(Ds, dim(left_virtualspace(ψ, 1))) + push!(E₀s, real(expectation_value(ψ, H, envs))) + push!(ΔEs, real(Es_step[1])) +end +[Ds E₀s ΔEs] +```` + +```` +9×3 Matrix{Float64}: + 165.0 -2.57351 -0.0418567 + 289.0 -2.58456 0.00418946 + 380.0 -2.58558 0.00073655 + 524.0 -2.58673 0.00949052 + 674.0 -2.58715 0.0111844 + 824.0 -2.58736 0.0116534 + 942.0 -2.58746 0.0122323 + 1058.0 -2.58757 0.0125425 + 1195.0 -2.58769 0.0128488 +```` + +````julia +fig = Figure() +ax = Axis(fig[1, 1]; xlabel = "D", ylabel = "ΔE / J", title = "[2 1 0], k = 2π/3") +scatterlines!(ax, Ds, ΔEs) +hlines!(ax, [0.0]; color = :gray, linestyle = :dash) +hlines!(ax, [0.0263]; color = :gray, linestyle = :dot) +fig +```` + +![](figure-1.png) + +As the bond dimension grows, the ground-state energy decreases, and the excitation energy turns positive already after the first expansion. +Apart from a dip at the second expansion, where VUMPS does not fully converge within `maxiter`, it then increases steadily. +It is still increasing at the largest bond dimension used here, and remains below the value ``ΔE / J = 0.0263`` reported in [devos2022](@cite) (dotted line). +Resolving this small gap therefore requires an even larger bond dimension. + +Finally, we compute the dispersion relation on top of the largest ground state, scanning half of the Brillouin zone. + +````julia +kspace = range(0, π, 10) +Es, _ = excitations(H, QuasiparticleAnsatz(), kspace, ψ, envs; sector, verbosity = 0) + +fig = Figure() +ax = Axis(fig[1, 1]; xlabel = "momentum", ylabel = "ΔE / J", title = "[2 1 0]") +lines!(ax, kspace, vec(real.(Es))) +vlines!(ax, [2π / 3]; color = :gray, linestyle = :dot) +fig +```` + +![](figure-2.png) + +The dispersion has its minimum at ``k = 2π/3``, where [devos2022](@cite) locates the gap of this model. + +--- + +*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* + diff --git a/docs/src/examples/excitations/su3-heisenberg/main.ipynb b/docs/src/examples/excitations/su3-heisenberg/main.ipynb new file mode 100644 index 000000000..9164610a4 --- /dev/null +++ b/docs/src/examples/excitations/su3-heisenberg/main.ipynb @@ -0,0 +1,306 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# The SU(3) Heisenberg chain\n", + "\n", + "The spin-1 Heisenberg chain has a gapped, symmetry-protected topological ground state, the Haldane phase.\n", + "In this example we study what happens when the SU(2) symmetry of the spins is promoted to SU(3), following [devos2022](@cite).\n", + "Every site of the chain carries the fully symmetric $[3 0 0]$ irrep of SU(3), which is ten-dimensional, and neighbouring sites are coupled by the SU(3)-invariant Heisenberg interaction\n", + "\n", + "$$\n", + "H = J ∑_{⟨i,j⟩} ∑_{a=1}^{8} Tᵃᵢ Tᵃⱼ\n", + "$$\n", + "\n", + "where the $Tᵃ$ are the eight generators of SU(3).\n", + "This is the direct analogue of the spin-1 chain, with the three generators of SU(2) replaced by the eight generators of SU(3).\n", + "\n", + "Just as in the Haldane gap example, imposing the full non-abelian symmetry block-diagonalizes every tensor, so that a modest number of variational parameters captures a large bond dimension.\n", + "The SU(3) Clebsch–Gordan coefficients are provided by [SUNRepresentations.jl](https://github.com/QuantumKitHub/SUNRepresentations.jl), whose `SU3Irrep` sectors plug directly into TensorKit.\n", + "\n", + "> **First-run cost**\n", + ">\n", + "> SUNRepresentations.jl computes Clebsch–Gordan coefficients on demand, and caches them both in memory and on disk.\n", + "> The first time a coefficient is needed, this computation can be slow.\n", + "\n", + "To follow this example you need the following packages:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using MPSKit, TensorKit, SUNRepresentations, CairoMakie" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For reproducibility of this page, we fix the seed of the random number generator:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using Random\n", + "Random.seed!(4321);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## The model\n", + "\n", + "Each site carries the $[3 0 0]$ irrep, which we specify through its highest weight." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "I300 = SU3Irrep(3, 0, 0)\n", + "Vphys = Vect[SU3Irrep](I300 => 1)\n", + "dim(I300)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To build the interaction, we use that it is diagonal in the total charge $c$ of a bond.\n", + "Writing $C₂ = ∑ₐ TᵃTᵃ$ for the quadratic Casimir, we have\n", + "\n", + "$$\n", + "∑ₐ Tᵃᵢ Tᵃⱼ = ½ (C₂(c) - 2 C₂([3 0 0])).\n", + "$$\n", + "\n", + "SUNRepresentations.jl provides the quadratic Casimir as `casimir(2, irrep)`, normalized such that the fundamental $[1 0 0]$ has $C₂ = 4/3$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "casimir(2, SU3Irrep(1, 0, 0)), casimir(2, I300)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The nearest-neighbour term is then a single operator that acts as a scalar on each fusion channel of $[3 0 0] ⊗ [3 0 0]$.\n", + "We construct it by filling in the blocks, one for every total charge $c$:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "function su3_heisenberg(Vphys; J = 1.0)\n", + " c2_site = casimir(2, SU3Irrep(3, 0, 0))\n", + " SS = zeros(ComplexF64, Vphys ⊗ Vphys ← Vphys ⊗ Vphys)\n", + " for (c, data) in blocks(SS)\n", + " coupling = (casimir(2, c) - 2 * c2_site) / 2\n", + " for i in axes(data, 1)\n", + " data[i, i] = J * coupling\n", + " end\n", + " end\n", + " return InfiniteMPOHamiltonian(SS)\n", + "end\n", + "\n", + "H = su3_heisenberg(Vphys)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The product $[3 0 0] ⊗ [3 0 0]$ decomposes into four irreps, and the antiferromagnetic coupling favours the channels with the lowest Casimir:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for (c, n) in directproduct(I300, I300)\n", + " coupling = (casimir(2, c) - 2 * casimir(2, I300)) / 2\n", + " println(\"$c:\\t dim = $(dim(c)),\\t coupling = $coupling\")\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Ground state\n", + "\n", + "The $[3 0 0]$ irrep has trivial triality, so that, as for the spin-1 chain, a translation-invariant state with a single-site unit cell is allowed.\n", + "The physical charge then cannot change the triality of the virtual spaces, so all virtual charges must have the same triality.\n", + "Here we use virtual charges with trivial triality, which give the lowest variational energy.\n", + "We start from a small virtual space, and optimize a random initial state with VUMPS:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Vvirt = Vect[SU3Irrep](\n", + " SU3Irrep(0, 0, 0) => 3,\n", + " SU3Irrep(2, 1, 0) => 6,\n", + " SU3Irrep(3, 0, 0) => 3,\n", + " SU3Irrep(3, 3, 0) => 3,\n", + " SU3Irrep(4, 2, 0) => 2,\n", + ")\n", + "ψ₀ = InfiniteMPS([Vphys], [Vvirt])\n", + "vumps = VUMPS(; maxiter = 400, tol = 1.0e-10, verbosity = 0)\n", + "ψ, envs, info = find_groundstate(ψ₀, H, vumps)\n", + "E₀ = real(expectation_value(ψ, H, envs))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Excitations\n", + "\n", + "On top of the ground state, we compute excitations with the quasiparticle ansatz.\n", + "Following [devos2022](@cite), we look for the lowest excitations in the adjoint $[2 1 0]$ sector, which has its minimum at momentum $k = 2π/3$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sector = SU3Irrep(2, 1, 0)\n", + "Es, _ = excitations(H, QuasiparticleAnsatz(), 2π / 3, ψ, envs; sector)\n", + "ΔE = real(Es[1])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The excitation energy comes out negative, which would mean that there are states below the ground state.\n", + "\n", + "> **Excitation energies are not variational**\n", + ">\n", + "> The quasiparticle ansatz returns the energy of the excited state *relative to the ground-state energy*.\n", + "> The ground state itself is only approximated, and its energy $E₀$ lies above the exact one.\n", + "> Any error in $E₀$ therefore directly shifts the excitation energies, and they are not an upper bound on the true gap.\n", + "> In particular, a gap that is small compared to the error in $E₀$ can come out too small, or even negative.\n", + "\n", + "The remedy is to improve the ground state.\n", + "Instead of guessing a larger virtual space, we let `changebonds` with `OptimalExpand` grow the bond dimension, by adding the directions that contribute the most to a two-site update of the converged state.\n", + "Here `truncrank` counts the total dimension that is added to each bond, including the dimension of the SU(3) irreps.\n", + "Every expanded state is then reconverged with VUMPS, and we recompute the excitation energy at $k = 2π/3$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "expand = OptimalExpand(; trunc = truncrank(150))\n", + "Ds = [dim(left_virtualspace(ψ, 1))]\n", + "E₀s = [E₀]\n", + "ΔEs = [ΔE]\n", + "for step in 1:8\n", + " global ψ, envs = changebonds(ψ, H, expand, envs)\n", + " global ψ, envs, info = find_groundstate(ψ, H, vumps, envs)\n", + " Es_step, _ = excitations(H, QuasiparticleAnsatz(), 2π / 3, ψ, envs; sector)\n", + " push!(Ds, dim(left_virtualspace(ψ, 1)))\n", + " push!(E₀s, real(expectation_value(ψ, H, envs)))\n", + " push!(ΔEs, real(Es_step[1]))\n", + "end\n", + "[Ds E₀s ΔEs]\n", + "fig = Figure()\n", + "ax = Axis(fig[1, 1]; xlabel = \"D\", ylabel = \"ΔE / J\", title = \"[2 1 0], k = 2π/3\")\n", + "scatterlines!(ax, Ds, ΔEs)\n", + "hlines!(ax, [0.0]; color = :gray, linestyle = :dash)\n", + "hlines!(ax, [0.0263]; color = :gray, linestyle = :dot)\n", + "fig" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As the bond dimension grows, the ground-state energy decreases, and the excitation energy turns positive already after the first expansion.\n", + "Apart from a dip at the second expansion, where VUMPS does not fully converge within `maxiter`, it then increases steadily.\n", + "It is still increasing at the largest bond dimension used here, and remains below the value $ΔE / J = 0.0263$ reported in [devos2022](@cite) (dotted line).\n", + "Resolving this small gap therefore requires an even larger bond dimension.\n", + "\n", + "Finally, we compute the dispersion relation on top of the largest ground state, scanning half of the Brillouin zone." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "kspace = range(0, π, 10)\n", + "Es, _ = excitations(H, QuasiparticleAnsatz(), kspace, ψ, envs; sector, verbosity = 0)\n", + "\n", + "fig = Figure()\n", + "ax = Axis(fig[1, 1]; xlabel = \"momentum\", ylabel = \"ΔE / J\", title = \"[2 1 0]\")\n", + "lines!(ax, kspace, vec(real.(Es)))\n", + "vlines!(ax, [2π / 3]; color = :gray, linestyle = :dot)\n", + "fig" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The dispersion has its minimum at $k = 2π/3$, where [devos2022](@cite) locates the gap of this model." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "\n", + "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.13.3", + "language": "julia", + "name": "julia-1.13" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.13.3" + } + }, + "nbformat": 4, + "nbformat_minor": 3 +} \ No newline at end of file diff --git a/docs/src/examples/groundstates/tfim-groundstate/figure-1.png b/docs/src/examples/groundstates/tfim-groundstate/figure-1.png new file mode 100644 index 000000000..0b9ba1f2b Binary files /dev/null and b/docs/src/examples/groundstates/tfim-groundstate/figure-1.png differ diff --git a/docs/src/examples/groundstates/tfim-groundstate/figure-2.png b/docs/src/examples/groundstates/tfim-groundstate/figure-2.png new file mode 100644 index 000000000..6ce16b5cd Binary files /dev/null and b/docs/src/examples/groundstates/tfim-groundstate/figure-2.png differ diff --git a/docs/src/examples/groundstates/tfim-groundstate/figure-3.png b/docs/src/examples/groundstates/tfim-groundstate/figure-3.png new file mode 100644 index 000000000..66d8aef2e Binary files /dev/null and b/docs/src/examples/groundstates/tfim-groundstate/figure-3.png differ diff --git a/docs/src/examples/groundstates/tfim-groundstate/index.md b/docs/src/examples/groundstates/tfim-groundstate/index.md new file mode 100644 index 000000000..6c463a1f9 --- /dev/null +++ b/docs/src/examples/groundstates/tfim-groundstate/index.md @@ -0,0 +1,184 @@ +```@meta +EditURL = "../../../../../examples/groundstates/tfim-groundstate/main.jl" +``` + +[![](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/QuantumKitHub/MPSKit.jl/gh-pages?filepath=dev/examples/groundstates/tfim-groundstate/main.ipynb) +[![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](https://nbviewer.jupyter.org/github/QuantumKitHub/MPSKit.jl/blob/gh-pages/dev/examples/groundstates/tfim-groundstate/main.ipynb) +[![](https://img.shields.io/badge/download-project-orange)](https://minhaskamal.github.io/DownGit/#/home?url=https://github.com/QuantumKitHub/MPSKit.jl/examples/tree/gh-pages/dev/examples/groundstates/tfim-groundstate) + +# The transverse-field Ising model + +In this example we locate the quantum phase transition of the transverse-field Ising model, by scanning the field across it and comparing a finite chain with a calculation directly in the thermodynamic limit. +The Hamiltonian is given by + +```math +H = -J (∑_{⟨i,j⟩} σᶻᵢ σᶻⱼ + g ∑ᵢ σˣᵢ) +``` + +where ``J`` sets the energy scale and the field ``g`` tunes the competition between the ``σᶻσᶻ`` interaction and the transverse ``σˣ`` field. +The model has a quantum critical point at ``g = 1``, which we will diagnose in three independent ways: the order parameter, the entanglement entropy and the correlation length. +To follow this example you need the following packages: + +````julia +using MPSKit, TensorKit, CairoMakie +using .ExampleModels +using TensorKitTensors.SpinOperators: σᶻ +```` + +For reproducibility of this page, we fix the seed of the random number generator: + +````julia +using Random +Random.seed!(1234); +```` + +We fix a finite chain length `L`, a bond dimension `D` and the set of field values to scan. +`D` is kept modest so the example runs quickly; increasing it sharpens the results in the thermodynamic limit. + +````julia +L = 16 +D = 8 +g_values = 0.1:0.1:2.0 +```` + +```` +0.1:0.1:2.0 +```` + +## Order parameter + +We compute the order parameter ``|⟨σᶻ⟩|`` at every field value, both on a finite chain and in the thermodynamic limit. + +For the finite chain we use an open chain of `L` sites and optimize with DMRG. +We average ``⟨σᶻᵢ⟩`` over the sites and take the absolute value, since a symmetry-broken state can have either sign. + +````julia +ψ₀_finite = FiniteMPS(L, ℂ^2, ℂ^D) +M_finite = map(g_values) do g + H = transverse_field_ising(; L, g) + ψ, envs, info = find_groundstate(ψ₀_finite, H, DMRG(; verbosity = 0)) + return abs(sum(i -> expectation_value(ψ, i => σᶻ()), 1:L)) / L +end +```` + +```` +20-element Vector{Float64}: + 0.9980617138856797 + 0.009304206174155815 + 9.886926889666883e-7 + 1.4949130880990819e-7 + 4.847831921565129e-9 + 2.4238082826113185e-10 + 5.090603802018776e-10 + 6.399930826746388e-12 + 1.1608529523047178e-12 + 9.505170272796344e-15 + 2.2076248830477028e-14 + 3.5562341351713215e-14 + 4.552007274692689e-15 + 2.311668137737342e-15 + 1.541625133765297e-15 + 6.59959564548599e-16 + 2.5681839634229746e-16 + 1.2731436158915082e-15 + 9.96912572737997e-16 + 1.9572671223156964e-16 +```` + +In the thermodynamic limit, we drop the length `L` to build the Hamiltonian on the infinite chain, and optimize an `InfiniteMPS` with VUMPS. +We keep the optimized states, since we will reuse them for the entropy and the correlation length below. +On the ordered side, the infinite state settles into one of the two symmetry-broken ground states, so we again take the absolute value. + +````julia +ψ₀_infinite = InfiniteMPS(ℂ^2, ℂ^D) +states_infinite = map(g_values) do g + H = transverse_field_ising(; g) + ψ, envs, info = find_groundstate(ψ₀_infinite, H, VUMPS(; verbosity = 0)) + return ψ +end +M_infinite = [abs(expectation_value(ψ, 1 => σᶻ())) for ψ in states_infinite] +```` + +```` +20-element Vector{Float64}: + 0.9987444968227268 + 0.9949102475942 + 0.9882803811214428 + 0.9784416036423226 + 0.9646786299603123 + 0.9457416090034619 + 0.9192767410625147 + 0.8801117382538193 + 0.8125392063206265 + 0.4785381092362469 + 1.4673123827080297e-12 + 9.027500968983304e-15 + 1.0614772949502084e-14 + 2.0311877180212447e-14 + 6.019490461639521e-15 + 7.421147030228778e-15 + 4.833806965809376e-15 + 4.174612044938186e-15 + 1.2895934320411578e-14 + 2.584737979205443e-16 +```` + +````julia +fig = Figure() +ax = Axis(fig[1, 1]; xlabel = "g", ylabel = "|⟨σᶻ⟩|") +scatter!(ax, g_values, M_finite; label = "L = $L") +scatter!(ax, g_values, M_infinite; label = "L = ∞") +vlines!(ax, [1.0]; color = :gray, linestyle = :dash) +axislegend(ax; position = :rt) +fig +```` + +![](figure-1.png) + +The two calculations agree in the disordered phase, but not in the ordered one. +The infinite state stays ordered up to ``g = 1`` and then collapses, locating the critical point. +On the finite chain, the magnetization vanishes well before ``g = 1``. +The exact ground state of a finite chain is the symmetric superposition of the two ordered states, which only needs a small bond dimension, so DMRG readily finds it. +Spontaneous symmetry breaking only happens in the thermodynamic limit, which is why the order parameter is best measured there. + +## Entanglement entropy + +At a critical point correlations become long-ranged, and the ground state is more entangled than deep in either phase. +For an `InfiniteMPS`, `entropy` returns the entanglement entropy of every bond in the unit cell, which here contains a single site. + +````julia +S_infinite = [real(only(entropy(ψ))) for ψ in states_infinite] + +fig, ax, _ = scatter(g_values, S_infinite; axis = (; xlabel = "g", ylabel = "S")) +vlines!(ax, [1.0]; color = :gray, linestyle = :dash) +fig +```` + +![](figure-2.png) + +The entropy peaks near ``g = 1``, the entanglement signature of the phase transition. + +## Correlation length + +Finally, the correlation length is extracted from the transfer matrix spectrum of the infinite state, and has no finite-chain analogue. +`correlation_length` returns one value per symmetry sector, and without symmetries there is only the trivial one. +It grows towards the critical point, so we plot it on a logarithmic scale. + +````julia +ξ_infinite = [only(values(correlation_length(ψ))) for ψ in states_infinite] + +fig, ax, _ = scatter(g_values, ξ_infinite; axis = (; xlabel = "g", ylabel = "ξ", yscale = log10)) +vlines!(ax, [1.0]; color = :gray, linestyle = :dash) +fig +```` + +![](figure-3.png) + +The correlation length peaks near ``g = 1`` as well. +At the critical point it should diverge, but an MPS with finite bond dimension `D` can only capture correlations up to a finite length, so the peak grows as `D` is increased. +The critical point itself is studied in more detail in the [Ising CFT](@ref "The Ising CFT spectrum") example, which matches its excitation spectrum to the predictions of conformal field theory. + +--- + +*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* + diff --git a/docs/src/examples/groundstates/tfim-groundstate/main.ipynb b/docs/src/examples/groundstates/tfim-groundstate/main.ipynb new file mode 100644 index 000000000..ccaa5d32f --- /dev/null +++ b/docs/src/examples/groundstates/tfim-groundstate/main.ipynb @@ -0,0 +1,329 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# The transverse-field Ising model\n", + "\n", + "In this example we locate the quantum phase transition of the transverse-field Ising model, by scanning the field across it and comparing a finite chain with a calculation directly in the thermodynamic limit.\n", + "The Hamiltonian is given by\n", + "\n", + "$$\n", + "H = -J (∑_{⟨i,j⟩} σᶻᵢ σᶻⱼ + g ∑ᵢ σˣᵢ)\n", + "$$\n", + "\n", + "where $J$ sets the energy scale and the field $g$ tunes the competition between the $σᶻσᶻ$ interaction and the transverse $σˣ$ field.\n", + "The model has a quantum critical point at $g = 1$, which we will diagnose in three independent ways: the order parameter, the entanglement entropy and the correlation length.\n", + "To follow this example you need the following packages:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using MPSKit, TensorKit, CairoMakie\n", + "# Model constructors for the MPSKit documentation and examples.\n", + "# Adapted from MPSKitModels.jl v0.4.7, src/models/{hamiltonians,transfermatrices}.jl:\n", + "# https://github.com/QuantumKitHub/MPSKitModels.jl/tree/v0.4.7/src/models\n", + "# Operators come from TensorKitTensors; chain assembly uses MPSKit directly.\n", + "#\n", + "# MIT License\n", + "# Copyright (c) 2021 Maarten Van Damme\n", + "#\n", + "# Permission is hereby granted, free of charge, to any person obtaining a copy\n", + "# of this software and associated documentation files (the \"Software\"), to deal\n", + "# in the Software without restriction, including without limitation the rights\n", + "# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell\n", + "# copies of the Software, and to permit persons to whom the Software is\n", + "# furnished to do so, subject to the following conditions:\n", + "#\n", + "# The above copyright notice and this permission notice shall be included in all\n", + "# copies or substantial portions of the Software.\n", + "#\n", + "# THE SOFTWARE IS PROVIDED \"AS IS\", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR\n", + "# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,\n", + "# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE\n", + "# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER\n", + "# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,\n", + "# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE\n", + "# SOFTWARE.\n", + "\n", + "module ExampleModels\n", + "\n", + "using MPSKit, TensorKit\n", + "using TensorKitTensors: SpinOperators, BosonOperators, HubbardOperators\n", + "\n", + "export chain_hamiltonian, transverse_field_ising, heisenberg_XXX\n", + "export hubbard_model, bose_hubbard_model, hard_hexagon\n", + "\n", + "# L specifies an open finite chain; otherwise unitcell specifies an infinite chain.\n", + "function chain_hamiltonian(twosite, onesite = nothing; L = nothing, unitcell::Integer = 1)\n", + " nsites = isnothing(L) ? unitcell : L\n", + " nsites isa Integer && nsites > 0 || throw(ArgumentError(\"chain length must be a positive integer\"))\n", + " nbonds = isnothing(L) ? nsites : nsites - 1\n", + " terms = Pair[(i, i + 1) => twosite for i in 1:nbonds]\n", + " if !isnothing(onesite)\n", + " append!(terms, [i => onesite for i in 1:nsites])\n", + " end\n", + " spaces = fill(space(twosite, 1), nsites)\n", + " constructor = isnothing(L) ? InfiniteMPOHamiltonian : FiniteMPOHamiltonian\n", + " return constructor(spaces, terms)\n", + "end\n", + "\n", + "# H = -J sum(σz_i σz_{i+1} + g σx_i).\n", + "function transverse_field_ising(\n", + " T::Type{<:Number} = ComplexF64, symmetry::Type{<:Sector} = Trivial;\n", + " J = 1.0, g = 1.0, kwargs...\n", + " )\n", + " ZZ = 4 * SpinOperators.S_z_S_z(T, symmetry; spin = 1 // 2)\n", + " X = SpinOperators.σˣ(T, symmetry)\n", + " return chain_hamiltonian(-J * ZZ, -J * g * X; kwargs...)\n", + "end\n", + "transverse_field_ising(symmetry::Type{<:Sector}; kwargs...) =\n", + " transverse_field_ising(ComplexF64, symmetry; kwargs...)\n", + "\n", + "# H = J sum(S_i ⋅ S_{i+1}); preserve the examples' spin-1 default.\n", + "function heisenberg_XXX(\n", + " T::Type{<:Number} = ComplexF64, symmetry::Type{<:Sector} = Trivial;\n", + " J = 1.0, spin = 1, kwargs...\n", + " )\n", + " term = J * SpinOperators.S_exchange(T, symmetry; spin)\n", + " return chain_hamiltonian(term; kwargs...)\n", + "end\n", + "heisenberg_XXX(symmetry::Type{<:Sector}; kwargs...) =\n", + " heisenberg_XXX(ComplexF64, symmetry; kwargs...)\n", + "\n", + "# H = -t sum(e⁺_i e⁻_{i+1} + h.c.) + U sum(n↑ n↓) - mu sum(n).\n", + "function hubbard_model(\n", + " T::Type{<:Number} = ComplexF64, particle_symmetry::Type{<:Sector} = Trivial,\n", + " spin_symmetry::Type{<:Sector} = Trivial;\n", + " t = 1.0, U = 1.0, mu = 0.0, kwargs...\n", + " )\n", + " # e_hopping includes the fermionic reordering sign and is Hermitian.\n", + " hopping = HubbardOperators.e_hopping(T, particle_symmetry, spin_symmetry)\n", + " interaction = HubbardOperators.ud_num(T, particle_symmetry, spin_symmetry)\n", + " N = HubbardOperators.e_num(T, particle_symmetry, spin_symmetry)\n", + " return chain_hamiltonian(-t * hopping, U * interaction - mu * N; kwargs...)\n", + "end\n", + "\n", + "# H = -t sum(b⁺_i b⁻_{i+1} + h.c.) + U/2 sum(n(n-1)) - mu sum(n).\n", + "function bose_hubbard_model(\n", + " T::Type{<:Number} = ComplexF64, symmetry::Type{<:Sector} = Trivial;\n", + " cutoff::Integer = 5, t = 1.0, U = 1.0, mu = 0.0, kwargs...\n", + " )\n", + " hopping = BosonOperators.b_plus_b_min(T, symmetry; cutoff) +\n", + " BosonOperators.b_min_b_plus(T, symmetry; cutoff)\n", + " N = BosonOperators.b_num(T, symmetry; cutoff)\n", + " interaction = N * (N - id(domain(N)))\n", + " return chain_hamiltonian(-t * hopping, U / 2 * interaction - mu * N; kwargs...)\n", + "end\n", + "\n", + "function hard_hexagon(T::Type{<:Number} = ComplexF64)\n", + " P = Vect[FibonacciAnyon](:τ => 1)\n", + " O = ones(T, P ⊗ P ← P ⊗ P)\n", + " block(O, FibonacciAnyon(:I)) .*= 0\n", + " return InfiniteMPO([O])\n", + "end\n", + "\n", + "end\n", + "\n", + "using .ExampleModels\n", + "using TensorKitTensors.SpinOperators: σᶻ" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For reproducibility of this page, we fix the seed of the random number generator:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using Random\n", + "Random.seed!(1234);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We fix a finite chain length `L`, a bond dimension `D` and the set of field values to scan.\n", + "`D` is kept modest so the example runs quickly; increasing it sharpens the results in the thermodynamic limit." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "L = 16\n", + "D = 8\n", + "g_values = 0.1:0.1:2.0" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Order parameter\n", + "\n", + "We compute the order parameter $|⟨σᶻ⟩|$ at every field value, both on a finite chain and in the thermodynamic limit.\n", + "\n", + "For the finite chain we use an open chain of `L` sites and optimize with DMRG.\n", + "We average $⟨σᶻᵢ⟩$ over the sites and take the absolute value, since a symmetry-broken state can have either sign." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ψ₀_finite = FiniteMPS(L, ℂ^2, ℂ^D)\n", + "M_finite = map(g_values) do g\n", + " H = transverse_field_ising(; L, g)\n", + " ψ, envs, info = find_groundstate(ψ₀_finite, H, DMRG(; verbosity = 0))\n", + " return abs(sum(i -> expectation_value(ψ, i => σᶻ()), 1:L)) / L\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the thermodynamic limit, we drop the length `L` to build the Hamiltonian on the infinite chain, and optimize an `InfiniteMPS` with VUMPS.\n", + "We keep the optimized states, since we will reuse them for the entropy and the correlation length below.\n", + "On the ordered side, the infinite state settles into one of the two symmetry-broken ground states, so we again take the absolute value." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ψ₀_infinite = InfiniteMPS(ℂ^2, ℂ^D)\n", + "states_infinite = map(g_values) do g\n", + " H = transverse_field_ising(; g)\n", + " ψ, envs, info = find_groundstate(ψ₀_infinite, H, VUMPS(; verbosity = 0))\n", + " return ψ\n", + "end\n", + "M_infinite = [abs(expectation_value(ψ, 1 => σᶻ())) for ψ in states_infinite]\n", + "fig = Figure()\n", + "ax = Axis(fig[1, 1]; xlabel = \"g\", ylabel = \"|⟨σᶻ⟩|\")\n", + "scatter!(ax, g_values, M_finite; label = \"L = $L\")\n", + "scatter!(ax, g_values, M_infinite; label = \"L = ∞\")\n", + "vlines!(ax, [1.0]; color = :gray, linestyle = :dash)\n", + "axislegend(ax; position = :rt)\n", + "fig" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The two calculations agree in the disordered phase, but not in the ordered one.\n", + "The infinite state stays ordered up to $g = 1$ and then collapses, locating the critical point.\n", + "On the finite chain, the magnetization vanishes well before $g = 1$.\n", + "The exact ground state of a finite chain is the symmetric superposition of the two ordered states, which only needs a small bond dimension, so DMRG readily finds it.\n", + "Spontaneous symmetry breaking only happens in the thermodynamic limit, which is why the order parameter is best measured there." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Entanglement entropy\n", + "\n", + "At a critical point correlations become long-ranged, and the ground state is more entangled than deep in either phase.\n", + "For an `InfiniteMPS`, `entropy` returns the entanglement entropy of every bond in the unit cell, which here contains a single site." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "S_infinite = [real(only(entropy(ψ))) for ψ in states_infinite]\n", + "\n", + "fig, ax, _ = scatter(g_values, S_infinite; axis = (; xlabel = \"g\", ylabel = \"S\"))\n", + "vlines!(ax, [1.0]; color = :gray, linestyle = :dash)\n", + "fig" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The entropy peaks near $g = 1$, the entanglement signature of the phase transition." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Correlation length\n", + "\n", + "Finally, the correlation length is extracted from the transfer matrix spectrum of the infinite state, and has no finite-chain analogue.\n", + "`correlation_length` returns one value per symmetry sector, and without symmetries there is only the trivial one.\n", + "It grows towards the critical point, so we plot it on a logarithmic scale." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ξ_infinite = [only(values(correlation_length(ψ))) for ψ in states_infinite]\n", + "\n", + "fig, ax, _ = scatter(g_values, ξ_infinite; axis = (; xlabel = \"g\", ylabel = \"ξ\", yscale = log10))\n", + "vlines!(ax, [1.0]; color = :gray, linestyle = :dash)\n", + "fig" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The correlation length peaks near $g = 1$ as well.\n", + "At the critical point it should diverge, but an MPS with finite bond dimension `D` can only capture correlations up to a finite length, so the peak grows as `D` is increased.\n", + "The critical point itself is studied in more detail in the Ising CFT example, which matches its excitation spectrum to the predictions of conformal field theory." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "---\n", + "\n", + "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 1.13.3", + "language": "julia", + "name": "julia-1.13" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.13.3" + } + }, + "nbformat": 4, + "nbformat_minor": 3 +} \ No newline at end of file diff --git a/docs/src/examples/index.md b/docs/src/examples/index.md index 1f91825d8..4dd2f69b1 100644 --- a/docs/src/examples/index.md +++ b/docs/src/examples/index.md @@ -21,6 +21,14 @@ The sidebar lists the same pages in the same order. ## Ground states +### [The transverse-field Ising model](groundstates/tfim-groundstate/index.md) + +![](groundstates/tfim-groundstate/figure-1.png) + +Scans the transverse field across the critical point of the Ising chain and locates the transition three independent ways: the order parameter on a finite chain and directly in the thermodynamic limit, the entanglement entropy, and the correlation length. +Combines the basic finite and infinite ground-state tools into one case study, without any symmetries. +**Complexity: introductory.** + ### [The Ising CFT spectrum](groundstates/ising-cft/index.md) ![](groundstates/ising-cft/figure-2.png) @@ -74,6 +82,15 @@ Computes the Haldane gap of the spin-1 Heisenberg antiferromagnet in two complem Introduces SU(2)-symmetric tensors for both finite and infinite MPS. **Complexity: intermediate.** +### [The SU(3) Heisenberg chain](excitations/su3-heisenberg/index.md) + +![](excitations/su3-heisenberg/figure-1.png) + +Generalizes the spin-1 Heisenberg chain to a chain of ``[3 0 0]`` irreps of SU(3), building the symmetric coupling from quadratic Casimirs and optimizing a uniform ground state with VUMPS. +Growing the bond dimension with `changebonds` shows that quasiparticle excitation energies are not variational: the small gap at ``k = 2π/3`` comes out negative for a poorly converged ground state, and increases as the ground state improves. +Shows that the SU(2) workflow of the Haldane gap example carries over to SU(N) symmetries through SUNRepresentations.jl. +**Complexity: advanced.** + ## Dynamics & finite temperature ### [DQPT in the Ising model](dynamics/ising-dqpt/index.md) diff --git a/examples/Cache.toml b/examples/Cache.toml index 6f2cc32b2..296a00e2c 100644 --- a/examples/Cache.toml +++ b/examples/Cache.toml @@ -1,16 +1,18 @@ [excitations] -haldane = "2e79dd49111097b0541bb6039e2c4e23f615b61d59611da6fd633df229a20b5f" +haldane = "4a4576752d9a394f18177d4884fa00c60f6be4fb20084e0007b0a7e34457c225" +su3-heisenberg = "610c3fde3b5d5340b8352a96f8ef77b92e125fb2445edfab72f1ae65b4088615" [statmech] -hard-hexagon = "68269e430555ac77e0ed1d1c3b44e553ffa3bfb6a37112eb8c2aff350cd81e1c" +hard-hexagon = "55853c4136b6c0b6b91741dc28ff06cb4b53f3a808ef34d3b6abc8f93c91bfbd" [dynamics] -ising-dqpt = "cd6d0545172ada42735a7847d7da35be46bdf7a1b038366c7fa8046585269c92" -xy-finiteT = "11e794f7e8700da83bbeb75da8c697a47957c9469b830496974c7df442d854ab" +ising-dqpt = "bd3313835872764e28f77750433ed72b8097ee74a295afa9947a151bc80be6ef" +xy-finiteT = "2263b2b3fd3ddaaf8164b2d4866088028ca535f925f37d13a77e4c4965099d97" [groundstates] -bose-hubbard = "5411a7d968e9fb0071528f81946536046d12820f39dc358aaae1e1a3be9580bf" -haldane-spt = "4ad3562ff82691ac31e0851c30bf5c7975d7bbcb06494f93ba168084ef912d42" -xxz-heisenberg = "aad20bd24766dce3ee103e815723d1d5c15951a6f69557660525308a94b5e229" -ising-cft = "9715c2348b7bf056bf1742d4f0940faa3eb3d30857d439e027993fb8fbb98e0d" -hubbard = "66faf9d13d42abae30248d288f1e628c5cc5664e0ee69c6faf34b6b06d5ef75c" +tfim-groundstate = "c68c90200aa2a11df7a5a13ff996cae42837e88fc10c7d126fa9d6971c4fac24" +bose-hubbard = "d2c28ca050a2f0a7d039ae145052b1fbfbcd2114bfd2f33bd4067526918dedc0" +haldane-spt = "ecafb4eb37da88444ef8ffad7c5d237eb8f6ff7c99eef6b6ce01dc2169798c33" +xxz-heisenberg = "ceab3a8607ad3ad28d536309997f8db1fbf4b732ecea119ccfb314b9799f7f80" +ising-cft = "4881eac3246e7a1027f028a2b31f1389ddd2b84ad780d603df5f80f7f08b8a99" +hubbard = "0dcfd341636c9acdadbd024de7f36523b8c74c4465e66796c3b6ee053fc64607" diff --git a/examples/Project.toml b/examples/Project.toml index 262cac131..6a88e78f4 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -10,6 +10,7 @@ Optim = "429524aa-4258-5aef-a3af-852621145aeb" Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" +SUNRepresentations = "1a50b95c-7aac-476d-a9ce-2bfc675fc617" TensorKit = "07d1fe3e-3e46-537d-9eac-e9e13d0d4cec" TensorKitTensors = "41b62e7d-e9d1-4e23-942c-79a97adf954b" TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" diff --git a/examples/excitations/su3-heisenberg/main.jl b/examples/excitations/su3-heisenberg/main.jl new file mode 100644 index 000000000..f10066800 --- /dev/null +++ b/examples/excitations/su3-heisenberg/main.jl @@ -0,0 +1,173 @@ +md""" +# The SU(3) Heisenberg chain + +The spin-1 Heisenberg chain has a gapped, symmetry-protected topological ground state, the Haldane phase. +In this example we study what happens when the SU(2) symmetry of the spins is promoted to SU(3), following [devos2022](@cite). +Every site of the chain carries the fully symmetric ``[3 0 0]`` irrep of SU(3), which is ten-dimensional, and neighbouring sites are coupled by the SU(3)-invariant Heisenberg interaction + +```math +H = J ∑_{⟨i,j⟩} ∑_{a=1}^{8} Tᵃᵢ Tᵃⱼ +``` + +where the ``Tᵃ`` are the eight generators of SU(3). +This is the direct analogue of the spin-1 chain, with the three generators of SU(2) replaced by the eight generators of SU(3). + +Just as in the [Haldane gap](@ref "The Haldane gap") example, imposing the full non-abelian symmetry block-diagonalizes every tensor, so that a modest number of variational parameters captures a large bond dimension. +The SU(3) Clebsch–Gordan coefficients are provided by [SUNRepresentations.jl](https://github.com/QuantumKitHub/SUNRepresentations.jl), whose `SU3Irrep` sectors plug directly into TensorKit. + +!!! warning "First-run cost" + SUNRepresentations.jl computes Clebsch–Gordan coefficients on demand, and caches them both in memory and on disk. + The first time a coefficient is needed, this computation can be slow. + +To follow this example you need the following packages: +""" + +using MPSKit, TensorKit, SUNRepresentations, CairoMakie + +# For reproducibility of this page, we fix the seed of the random number generator: + +using Random +Random.seed!(4321); + +md""" +## The model + +Each site carries the ``[3 0 0]`` irrep, which we specify through its highest weight. +""" + +I300 = SU3Irrep(3, 0, 0) +Vphys = Vect[SU3Irrep](I300 => 1) +dim(I300) + +md""" +To build the interaction, we use that it is diagonal in the total charge ``c`` of a bond. +Writing ``C₂ = ∑ₐ TᵃTᵃ`` for the quadratic Casimir, we have + +```math +∑ₐ Tᵃᵢ Tᵃⱼ = ½ (C₂(c) - 2 C₂([3 0 0])). +``` + +SUNRepresentations.jl provides the quadratic Casimir as `casimir(2, irrep)`, normalized such that the fundamental ``[1 0 0]`` has ``C₂ = 4/3``. +""" + +casimir(2, SU3Irrep(1, 0, 0)), casimir(2, I300) + +md""" +The nearest-neighbour term is then a single operator that acts as a scalar on each fusion channel of ``[3 0 0] ⊗ [3 0 0]``. +We construct it by filling in the blocks, one for every total charge ``c``: +""" + +function su3_heisenberg(Vphys; J = 1.0) + c2_site = casimir(2, SU3Irrep(3, 0, 0)) + SS = zeros(ComplexF64, Vphys ⊗ Vphys ← Vphys ⊗ Vphys) + for (c, data) in blocks(SS) + coupling = (casimir(2, c) - 2 * c2_site) / 2 + for i in axes(data, 1) + data[i, i] = J * coupling + end + end + return InfiniteMPOHamiltonian(SS) +end + +H = su3_heisenberg(Vphys) + +md""" +The product ``[3 0 0] ⊗ [3 0 0]`` decomposes into four irreps, and the antiferromagnetic coupling favours the channels with the lowest Casimir: +""" + +for (c, n) in directproduct(I300, I300) + coupling = (casimir(2, c) - 2 * casimir(2, I300)) / 2 + println("$c:\t dim = $(dim(c)),\t coupling = $coupling") +end + +md""" +## Ground state + +The ``[3 0 0]`` irrep has trivial triality, so that, as for the spin-1 chain, a translation-invariant state with a single-site unit cell is allowed. +The physical charge then cannot change the triality of the virtual spaces, so all virtual charges must have the same triality. +Here we use virtual charges with trivial triality, which give the lowest variational energy. +We start from a small virtual space, and optimize a random initial state with VUMPS: +""" + +Vvirt = Vect[SU3Irrep]( + SU3Irrep(0, 0, 0) => 3, + SU3Irrep(2, 1, 0) => 6, + SU3Irrep(3, 0, 0) => 3, + SU3Irrep(3, 3, 0) => 3, + SU3Irrep(4, 2, 0) => 2, +) +ψ₀ = InfiniteMPS([Vphys], [Vvirt]) +vumps = VUMPS(; maxiter = 400, tol = 1.0e-10, verbosity = 0) +ψ, envs, info = find_groundstate(ψ₀, H, vumps) +E₀ = real(expectation_value(ψ, H, envs)) + +md""" +## Excitations + +On top of the ground state, we compute excitations with the quasiparticle ansatz. +Following [devos2022](@cite), we look for the lowest excitations in the adjoint ``[2 1 0]`` sector, which has its minimum at momentum ``k = 2π/3``. +""" + +sector = SU3Irrep(2, 1, 0) +Es, _ = excitations(H, QuasiparticleAnsatz(), 2π / 3, ψ, envs; sector) +ΔE = real(Es[1]) + +md""" +The excitation energy comes out negative, which would mean that there are states below the ground state. + +!!! warning "Excitation energies are not variational" + The quasiparticle ansatz returns the energy of the excited state *relative to the ground-state energy*. + The ground state itself is only approximated, and its energy ``E₀`` lies above the exact one. + Any error in ``E₀`` therefore directly shifts the excitation energies, and they are not an upper bound on the true gap. + In particular, a gap that is small compared to the error in ``E₀`` can come out too small, or even negative. + +The remedy is to improve the ground state. +Instead of guessing a larger virtual space, we let `changebonds` with `OptimalExpand` grow the bond dimension, by adding the directions that contribute the most to a two-site update of the converged state. +Here `truncrank` counts the total dimension that is added to each bond, including the dimension of the SU(3) irreps. +Every expanded state is then reconverged with VUMPS, and we recompute the excitation energy at ``k = 2π/3``. +""" + +expand = OptimalExpand(; trunc = truncrank(150)) +Ds = [dim(left_virtualspace(ψ, 1))] +E₀s = [E₀] +ΔEs = [ΔE] +for step in 1:8 + global ψ, envs = changebonds(ψ, H, expand, envs) + global ψ, envs, info = find_groundstate(ψ, H, vumps, envs) + Es_step, _ = excitations(H, QuasiparticleAnsatz(), 2π / 3, ψ, envs; sector) + push!(Ds, dim(left_virtualspace(ψ, 1))) + push!(E₀s, real(expectation_value(ψ, H, envs))) + push!(ΔEs, real(Es_step[1])) +end +[Ds E₀s ΔEs] + +#+ + +fig = Figure() +ax = Axis(fig[1, 1]; xlabel = "D", ylabel = "ΔE / J", title = "[2 1 0], k = 2π/3") +scatterlines!(ax, Ds, ΔEs) +hlines!(ax, [0.0]; color = :gray, linestyle = :dash) +hlines!(ax, [0.0263]; color = :gray, linestyle = :dot) +fig + +md""" +As the bond dimension grows, the ground-state energy decreases, and the excitation energy turns positive already after the first expansion. +Apart from a dip at the second expansion, where VUMPS does not fully converge within `maxiter`, it then increases steadily. +It is still increasing at the largest bond dimension used here, and remains below the value ``ΔE / J = 0.0263`` reported in [devos2022](@cite) (dotted line). +Resolving this small gap therefore requires an even larger bond dimension. + +Finally, we compute the dispersion relation on top of the largest ground state, scanning half of the Brillouin zone. +""" + +kspace = range(0, π, 10) +Es, _ = excitations(H, QuasiparticleAnsatz(), kspace, ψ, envs; sector, verbosity = 0) + +fig = Figure() +ax = Axis(fig[1, 1]; xlabel = "momentum", ylabel = "ΔE / J", title = "[2 1 0]") +lines!(ax, kspace, vec(real.(Es))) +vlines!(ax, [2π / 3]; color = :gray, linestyle = :dot) +fig + +md""" +The dispersion has its minimum at ``k = 2π/3``, where [devos2022](@cite) locates the gap of this model. +""" diff --git a/examples/groundstates/tfim-groundstate/main.jl b/examples/groundstates/tfim-groundstate/main.jl new file mode 100644 index 000000000..e14035d76 --- /dev/null +++ b/examples/groundstates/tfim-groundstate/main.jl @@ -0,0 +1,118 @@ +md""" +# The transverse-field Ising model + +In this example we locate the quantum phase transition of the transverse-field Ising model, by scanning the field across it and comparing a finite chain with a calculation directly in the thermodynamic limit. +The Hamiltonian is given by + +```math +H = -J (∑_{⟨i,j⟩} σᶻᵢ σᶻⱼ + g ∑ᵢ σˣᵢ) +``` + +where ``J`` sets the energy scale and the field ``g`` tunes the competition between the ``σᶻσᶻ`` interaction and the transverse ``σˣ`` field. +The model has a quantum critical point at ``g = 1``, which we will diagnose in three independent ways: the order parameter, the entanglement entropy and the correlation length. +To follow this example you need the following packages: +""" + +using MPSKit, TensorKit, CairoMakie +include(joinpath(@__DIR__, "..", "..", "..", "docs", "src", "assets", "models.jl")) #hide +using .ExampleModels +using TensorKitTensors.SpinOperators: σᶻ + +# For reproducibility of this page, we fix the seed of the random number generator: + +using Random +Random.seed!(1234); + +md""" +We fix a finite chain length `L`, a bond dimension `D` and the set of field values to scan. +`D` is kept modest so the example runs quickly; increasing it sharpens the results in the thermodynamic limit. +""" + +L = 16 +D = 8 +g_values = 0.1:0.1:2.0 + +md""" +## Order parameter + +We compute the order parameter ``|⟨σᶻ⟩|`` at every field value, both on a finite chain and in the thermodynamic limit. + +For the finite chain we use an open chain of `L` sites and optimize with DMRG. +We average ``⟨σᶻᵢ⟩`` over the sites and take the absolute value, since a symmetry-broken state can have either sign. +""" + +ψ₀_finite = FiniteMPS(L, ℂ^2, ℂ^D) +M_finite = map(g_values) do g + H = transverse_field_ising(; L, g) + ψ, envs, info = find_groundstate(ψ₀_finite, H, DMRG(; verbosity = 0)) + return abs(sum(i -> expectation_value(ψ, i => σᶻ()), 1:L)) / L +end + +md""" +In the thermodynamic limit, we drop the length `L` to build the Hamiltonian on the infinite chain, and optimize an `InfiniteMPS` with VUMPS. +We keep the optimized states, since we will reuse them for the entropy and the correlation length below. +On the ordered side, the infinite state settles into one of the two symmetry-broken ground states, so we again take the absolute value. +""" + +ψ₀_infinite = InfiniteMPS(ℂ^2, ℂ^D) +states_infinite = map(g_values) do g + H = transverse_field_ising(; g) + ψ, envs, info = find_groundstate(ψ₀_infinite, H, VUMPS(; verbosity = 0)) + return ψ +end +M_infinite = [abs(expectation_value(ψ, 1 => σᶻ())) for ψ in states_infinite] + +#+ + +fig = Figure() +ax = Axis(fig[1, 1]; xlabel = "g", ylabel = "|⟨σᶻ⟩|") +scatter!(ax, g_values, M_finite; label = "L = $L") +scatter!(ax, g_values, M_infinite; label = "L = ∞") +vlines!(ax, [1.0]; color = :gray, linestyle = :dash) +axislegend(ax; position = :rt) +fig + +md""" +The two calculations agree in the disordered phase, but not in the ordered one. +The infinite state stays ordered up to ``g = 1`` and then collapses, locating the critical point. +On the finite chain, the magnetization vanishes well before ``g = 1``. +The exact ground state of a finite chain is the symmetric superposition of the two ordered states, which only needs a small bond dimension, so DMRG readily finds it. +Spontaneous symmetry breaking only happens in the thermodynamic limit, which is why the order parameter is best measured there. +""" + +md""" +## Entanglement entropy + +At a critical point correlations become long-ranged, and the ground state is more entangled than deep in either phase. +For an `InfiniteMPS`, `entropy` returns the entanglement entropy of every bond in the unit cell, which here contains a single site. +""" + +S_infinite = [real(only(entropy(ψ))) for ψ in states_infinite] + +fig, ax, _ = scatter(g_values, S_infinite; axis = (; xlabel = "g", ylabel = "S")) +vlines!(ax, [1.0]; color = :gray, linestyle = :dash) +fig + +md""" +The entropy peaks near ``g = 1``, the entanglement signature of the phase transition. +""" + +md""" +## Correlation length + +Finally, the correlation length is extracted from the transfer matrix spectrum of the infinite state, and has no finite-chain analogue. +`correlation_length` returns one value per symmetry sector, and without symmetries there is only the trivial one. +It grows towards the critical point, so we plot it on a logarithmic scale. +""" + +ξ_infinite = [only(values(correlation_length(ψ))) for ψ in states_infinite] + +fig, ax, _ = scatter(g_values, ξ_infinite; axis = (; xlabel = "g", ylabel = "ξ", yscale = log10)) +vlines!(ax, [1.0]; color = :gray, linestyle = :dash) +fig + +md""" +The correlation length peaks near ``g = 1`` as well. +At the critical point it should diverge, but an MPS with finite bond dimension `D` can only capture correlations up to a finite length, so the peak grows as `D` is increased. +The critical point itself is studied in more detail in the [Ising CFT](@ref "The Ising CFT spectrum") example, which matches its excitation spectrum to the predictions of conformal field theory. +"""