diff --git a/docs/src/examples/excitations/1.su3-heisenberg/figure-1.png b/docs/src/examples/excitations/1.su3-heisenberg/figure-1.png new file mode 100644 index 000000000..ac9a559d5 Binary files /dev/null and b/docs/src/examples/excitations/1.su3-heisenberg/figure-1.png differ diff --git a/docs/src/examples/excitations/1.su3-heisenberg/index.md b/docs/src/examples/excitations/1.su3-heisenberg/index.md new file mode 100644 index 000000000..9b8bb1f40 --- /dev/null +++ b/docs/src/examples/excitations/1.su3-heisenberg/index.md @@ -0,0 +1,262 @@ +```@meta +EditURL = "../../../../../examples/excitations/1.su3-heisenberg/main.jl" +``` + +[![](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/QuantumKitHub/MPSKit.jl/gh-pages?filepath=dev/examples/excitations/1.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/1.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/1.su3-heisenberg) + +# The SU(3) Heisenberg chain + +The spin-1 Heisenberg chain famously has a gapped, symmetry-protected topological ground +state: the Haldane phase. +A natural question is what happens when the ``SU(2)`` symmetry of the spins is promoted to a +larger ``SU(N)`` symmetry. +In this example we study the ``SU(3)`` generalization introduced in [devos2022](@cite): a +uniform chain where every site carries the fully-symmetric ``[3\,0\,0]`` irrep of ``SU(3)`` +(the ten-dimensional representation), coupled by a nearest-neighbour ``SU(3)``-invariant +Heisenberg interaction + +```math +H = J \sum_{\langle i, j \rangle} \sum_{a=1}^{8} T^a_i T^a_j , +``` + +where the ``T^a`` are the eight generators of ``SU(3)``. +This is the direct analog of the spin-1 chain, with the three generators of ``SU(2)`` +replaced by the eight generators of ``SU(3)``. +The goal of this example is to exercise the full ``SU(N)`` toolchain end to end — building the +symmetric Hamiltonian, optimizing a uniform ground state with VUMPS, and scanning an +excitation dispersion — rather than to pin down the model's (very small) gap to high accuracy; +we return to that point in the closing discussion. + +Working with the full non-abelian ``SU(3)`` symmetry is what makes this tractable. +Just as ``SU(2)`` symmetry drastically shrinks the spin-1 problem (see +[Symmetries](@ref concept_symmetries)), imposing ``SU(3)`` block-diagonalizes every tensor +along its irreducible representations, so that a modest number of variational parameters +already captures a large physical bond dimension. +The Clebsch–Gordan coefficients that this requires are provided by +[SUNRepresentations.jl](https://github.com/QuantumKitHub/SUNRepresentations.jl), whose +`SU3Irrep` sectors plug directly into TensorKit's graded vector spaces. + +!!! warning "First-run cost" + SUNRepresentations.jl computes ``SU(N)`` Clebsch–Gordan coefficients on demand and + caches them both in memory and on disk (via `Scratch.jl`). + The **first** time a given coefficient is needed the computation can be slow; later runs + reuse the on-disk cache. + On a machine without a persisted scratch cache — such as a fresh continuous-integration + runner — this example therefore pays the full coefficient-generation cost once. + +````julia +using TensorKit +using SUNRepresentations +using MPSKit +using Plots +```` + +```` +Precompiling packages... + 960.3 ms ✓ ChunkCodecLibZstd + 1347.1 ms ✓ SparseArrayKit + 3750.1 ms ✓ FileIO + 25034.5 ms ✓ JLD2 + 7605.0 ms ✓ SUNRepresentations + 5 dependencies successfully precompiled in 38 seconds. 59 already precompiled. +Precompiling packages... + 2362.7 ms ✓ SparseArrayKit → SparseArrayKitSparseArrays + 1 dependency successfully precompiled in 3 seconds. 20 already precompiled. +Precompiling packages... + 1955.5 ms ✓ FileIO → HTTPExt + 4287.5 ms ✓ Plots → FileIOExt + 2 dependencies successfully precompiled in 8 seconds. 181 already precompiled. + +```` + +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 obtain from its highest weight and +represent as a one-dimensional graded space. + +````julia +I300 = SU3Irrep(3, 0, 0) +Vphys = Vect[SU3Irrep](I300 => 1) +@show dim(I300) +```` + +```` +10 +```` + +To build the interaction we use the same Casimir trick as for the ``SU(2)`` chain. +Writing ``C_2 = \sum_a T^a T^a`` for the quadratic Casimir operator, the two-site coupling is +diagonal in the total ``SU(3)`` charge ``c`` of a bond: + +```math +\sum_a T^a_i T^a_j = \tfrac{1}{2}\left( C_2(c) - C_2([3\,0\,0]) - C_2([3\,0\,0]) \right). +``` + +SUNRepresentations.jl provides the quadratic Casimir through `casimir(2, irrep)`. +Its normalization is the standard one with ``\mathrm{tr}(T^a T^b) = \tfrac{1}{2}\delta^{ab}``: +the fundamental ``[1\,0\,0]`` has ``C_2 = 4/3 = (N^2-1)/(2N)``. + +````julia +@show casimir(2, SU3Irrep(1, 0, 0)) +@show casimir(2, I300) +```` + +```` +6//1 +```` + +The nearest-neighbour term is then a single TensorKit operator that acts as a scalar on each +fusion channel of ``[3\,0\,0] \otimes [3\,0\,0]``. +Constructing it by iterating over the blocks — one block per total charge ``c`` — mirrors the +``SU(2)`` construction, but now with `casimir(2, ⋅)` supplying the eigenvalues. + +````julia +function su3_heisenberg(Vphys; J = 1.0) + site = SU3Irrep(3, 0, 0) + c2_site = casimir(2, site) + 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)) +| ⋮ + +```` + +We can read off the physics directly from the fusion channels. +The product ``[3\,0\,0] \otimes [3\,0\,0]`` decomposes into four irreps, and the +antiferromagnetic coupling (``J > 0``) favours the channels of lowest Casimir: + +````julia +for (c, n) in directproduct(I300, I300) + coupling = (casimir(2, c) - 2 * casimir(2, I300)) / 2 + println(rpad(string(c), 22), " dim = ", rpad(dim(c), 4), " coupling = ", coupling) +end +```` + +```` +Irrep[SU₃]((6, 0, 0)) dim = 28 coupling = 3//1 +Irrep[SU₃]((5, 1, 0)) dim = 35 coupling = 0//1 +Irrep[SU₃]((4, 2, 0)) dim = 27 coupling = -2//1 +Irrep[SU₃]((3, 3, 0)) dim = 10 coupling = -3//1 + +```` + +## Ground state + +The ``[3\,0\,0]`` irrep has trivial triality (congruency class ``0``), so — like the +integer-spin Haldane chain — a translation-invariant state on a single-site unit cell is +admissible. +Because the physical charge has triality ``0`` it cannot shift the triality of the virtual +bonds, so every virtual space must sit in a single congruency class; here we use the trivial +class ``0``, which we found to give the lowest variational energy. +Following [devos2022](@cite), we seed the optimization with a random ``SU(3)``-symmetric +uniform state and optimize with VUMPS. +The bond dimension is kept deliberately modest here so that the page builds quickly. + +````julia +Vvirt = Vect[SU3Irrep]( + SU3Irrep(0, 0, 0) => 8, + SU3Irrep(2, 1, 0) => 6, + SU3Irrep(3, 0, 0) => 3, + SU3Irrep(3, 3, 0) => 3, + SU3Irrep(4, 2, 0) => 2, +) +@show dim(Vvirt) + +ψ₀ = InfiniteMPS([Vphys], [Vvirt]) +ψ, envs, δ = find_groundstate(ψ₀, H, VUMPS(; maxiter = 200, tol = 1.0e-10, verbosity = 1)) +E₀ = real(expectation_value(ψ, H)) +println("ground-state energy per site: E₀ = $E₀") +```` + +```` +dim(Vvirt) = 170 +┌ Warning: Constructing an MPS from tensors that are not full rank +└ @ MPSKit src/states/infinitemps.jl:188 +ground-state energy per site: E₀ = -2.573505652431292 + +```` + +## The excitation spectrum + +On top of the ground state we compute excitations with the quasiparticle ansatz. +In the thermodynamic limit each excitation carries a definite momentum ``k`` and a definite +``SU(3)`` charge. +Following [devos2022](@cite), we look in the adjoint ``[2\,1\,0]`` sector, where that reference +finds the lowest branch, and scan across half of the Brillouin zone (the other half follows by +reflection). + +````julia +sector = SU3Irrep(2, 1, 0) +kspace = range(0, π, 10) +Es, _ = excitations(H, QuasiparticleAnsatz(), kspace, ψ, envs; sector, verbosity = 0) + +Δ, idx = findmin(real.(Es)) +println("branch minimum ΔE/J = $Δ at k = $(kspace[idx])") +```` + +```` +branch minimum ΔE/J = -0.041856660250395326 at k = 2.0943951023931953 + +```` + +````julia +plot( + kspace, real.(Es); + xaxis = "momentum k", yaxis = "ΔE / J", label = "[2 1 0] branch", + title = "SU(3) [3 0 0] excitation dispersion" +) +hline!([0.0]; color = :gray, linestyle = :dash, label = "") +vline!([2π / 3]; color = :gray, linestyle = :dot, label = "k = 2π/3") +```` + +![](figure-1.png) + +The branch has a pronounced soft mode at ``k = 2\pi/3``, exactly the momentum at which +[devos2022](@cite) locates the minimum of the dispersion, where that reference reports a gap +of ``\Delta / J = 0.0263``. + +The soft-mode *location* is reproduced cleanly, but the gap *value* is not, and it is worth +being precise about why. +The quasiparticle energies plotted above are a genuine variational upper bound on the +excitation energies only when they are measured relative to the *exact* ground state. +Here the ground state is a deliberately lightweight, modest-bond-dimension uniform ``SU(3)`` +MPS, and for this near-critical model it is not converged to that accuracy. +Because the reference state sits slightly too high in energy, the ansatz can find "excitations" +that lie below it, so near ``k = 2\pi/3`` the computed branch drops to essentially zero and even +slightly negative. +The dispersion here is therefore *not* a variational estimate of the true gap: the small value +``\Delta / J = 0.0263`` is quoted from [devos2022](@cite) and would require a carefully +converged, large-bond-dimension calculation to reproduce. +What this example does show is the complete ``SU(3)``-symmetric workflow — Hamiltonian +construction, VUMPS ground state, and a momentum-resolved excitation scan — with the soft mode +correctly located at ``k = 2\pi/3``. + +--- + +*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* + diff --git a/docs/src/examples/excitations/1.su3-heisenberg/main.ipynb b/docs/src/examples/excitations/1.su3-heisenberg/main.ipynb new file mode 100644 index 000000000..7c9958a7f --- /dev/null +++ b/docs/src/examples/excitations/1.su3-heisenberg/main.ipynb @@ -0,0 +1,302 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# The SU(3) Heisenberg chain\n", + "\n", + "The spin-1 Heisenberg chain famously has a gapped, symmetry-protected topological ground\n", + "state: the Haldane phase.\n", + "A natural question is what happens when the $SU(2)$ symmetry of the spins is promoted to a\n", + "larger $SU(N)$ symmetry.\n", + "In this example we study the $SU(3)$ generalization introduced in [devos2022](@cite): a\n", + "uniform chain where every site carries the fully-symmetric $[3\\,0\\,0]$ irrep of $SU(3)$\n", + "(the ten-dimensional representation), coupled by a nearest-neighbour $SU(3)$-invariant\n", + "Heisenberg interaction\n", + "\n", + "$$\n", + "H = J \\sum_{\\langle i, j \\rangle} \\sum_{a=1}^{8} T^a_i T^a_j ,\n", + "$$\n", + "\n", + "where the $T^a$ are the eight generators of $SU(3)$.\n", + "This is the direct analog of the spin-1 chain, with the three generators of $SU(2)$\n", + "replaced by the eight generators of $SU(3)$.\n", + "The goal of this example is to exercise the full $SU(N)$ toolchain end to end — building the\n", + "symmetric Hamiltonian, optimizing a uniform ground state with VUMPS, and scanning an\n", + "excitation dispersion — rather than to pin down the model's (very small) gap to high accuracy;\n", + "we return to that point in the closing discussion.\n", + "\n", + "Working with the full non-abelian $SU(3)$ symmetry is what makes this tractable.\n", + "Just as $SU(2)$ symmetry drastically shrinks the spin-1 problem (see\n", + "Symmetries), imposing $SU(3)$ block-diagonalizes every tensor\n", + "along its irreducible representations, so that a modest number of variational parameters\n", + "already captures a large physical bond dimension.\n", + "The Clebsch–Gordan coefficients that this requires are provided by\n", + "[SUNRepresentations.jl](https://github.com/QuantumKitHub/SUNRepresentations.jl), whose\n", + "`SU3Irrep` sectors plug directly into TensorKit's graded vector spaces.\n", + "\n", + "> **First-run cost**\n", + ">\n", + "> SUNRepresentations.jl computes $SU(N)$ Clebsch–Gordan coefficients on demand and\n", + "> caches them both in memory and on disk (via `Scratch.jl`).\n", + "> The **first** time a given coefficient is needed the computation can be slow; later runs\n", + "> reuse the on-disk cache.\n", + "> On a machine without a persisted scratch cache — such as a fresh continuous-integration\n", + "> runner — this example therefore pays the full coefficient-generation cost once." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using TensorKit\n", + "using SUNRepresentations\n", + "using MPSKit\n", + "using Plots" + ] + }, + { + "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 obtain from its highest weight and\n", + "represent as a one-dimensional graded space." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "I300 = SU3Irrep(3, 0, 0)\n", + "Vphys = Vect[SU3Irrep](I300 => 1)\n", + "@show dim(I300)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To build the interaction we use the same Casimir trick as for the $SU(2)$ chain.\n", + "Writing $C_2 = \\sum_a T^a T^a$ for the quadratic Casimir operator, the two-site coupling is\n", + "diagonal in the total $SU(3)$ charge $c$ of a bond:\n", + "\n", + "$$\n", + "\\sum_a T^a_i T^a_j = \\tfrac{1}{2}\\left( C_2(c) - C_2([3\\,0\\,0]) - C_2([3\\,0\\,0]) \\right).\n", + "$$\n", + "\n", + "SUNRepresentations.jl provides the quadratic Casimir through `casimir(2, irrep)`.\n", + "Its normalization is the standard one with $\\mathrm{tr}(T^a T^b) = \\tfrac{1}{2}\\delta^{ab}$:\n", + "the fundamental $[1\\,0\\,0]$ has $C_2 = 4/3 = (N^2-1)/(2N)$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "@show casimir(2, SU3Irrep(1, 0, 0))\n", + "@show casimir(2, I300)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The nearest-neighbour term is then a single TensorKit operator that acts as a scalar on each\n", + "fusion channel of $[3\\,0\\,0] \\otimes [3\\,0\\,0]$.\n", + "Constructing it by iterating over the blocks — one block per total charge $c$ — mirrors the\n", + "$SU(2)$ construction, but now with `casimir(2, ⋅)` supplying the eigenvalues." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "function su3_heisenberg(Vphys; J = 1.0)\n", + " site = SU3Irrep(3, 0, 0)\n", + " c2_site = casimir(2, site)\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": [ + "We can read off the physics directly from the fusion channels.\n", + "The product $[3\\,0\\,0] \\otimes [3\\,0\\,0]$ decomposes into four irreps, and the\n", + "antiferromagnetic coupling ($J > 0$) favours the channels of 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(rpad(string(c), 22), \" dim = \", rpad(dim(c), 4), \" coupling = \", coupling)\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Ground state\n", + "\n", + "The $[3\\,0\\,0]$ irrep has trivial triality (congruency class $0$), so — like the\n", + "integer-spin Haldane chain — a translation-invariant state on a single-site unit cell is\n", + "admissible.\n", + "Because the physical charge has triality $0$ it cannot shift the triality of the virtual\n", + "bonds, so every virtual space must sit in a single congruency class; here we use the trivial\n", + "class $0$, which we found to give the lowest variational energy.\n", + "Following [devos2022](@cite), we seed the optimization with a random $SU(3)$-symmetric\n", + "uniform state and optimize with VUMPS.\n", + "The bond dimension is kept deliberately modest here so that the page builds quickly." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Vvirt = Vect[SU3Irrep](\n", + " SU3Irrep(0, 0, 0) => 8,\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", + "@show dim(Vvirt)\n", + "\n", + "ψ₀ = InfiniteMPS([Vphys], [Vvirt])\n", + "ψ, envs, δ = find_groundstate(ψ₀, H, VUMPS(; maxiter = 200, tol = 1.0e-10, verbosity = 1))\n", + "E₀ = real(expectation_value(ψ, H))\n", + "println(\"ground-state energy per site: E₀ = $E₀\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## The excitation spectrum\n", + "\n", + "On top of the ground state we compute excitations with the quasiparticle ansatz.\n", + "In the thermodynamic limit each excitation carries a definite momentum $k$ and a definite\n", + "$SU(3)$ charge.\n", + "Following [devos2022](@cite), we look in the adjoint $[2\\,1\\,0]$ sector, where that reference\n", + "finds the lowest branch, and scan across half of the Brillouin zone (the other half follows by\n", + "reflection)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sector = SU3Irrep(2, 1, 0)\n", + "kspace = range(0, π, 10)\n", + "Es, _ = excitations(H, QuasiparticleAnsatz(), kspace, ψ, envs; sector, verbosity = 0)\n", + "\n", + "Δ, idx = findmin(real.(Es))\n", + "println(\"branch minimum ΔE/J = $Δ at k = $(kspace[idx])\")\n", + "plot(\n", + " kspace, real.(Es);\n", + " xaxis = \"momentum k\", yaxis = \"ΔE / J\", label = \"[2 1 0] branch\",\n", + " title = \"SU(3) [3 0 0] excitation dispersion\"\n", + ")\n", + "hline!([0.0]; color = :gray, linestyle = :dash, label = \"\")\n", + "vline!([2π / 3]; color = :gray, linestyle = :dot, label = \"k = 2π/3\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The branch has a pronounced soft mode at $k = 2\\pi/3$, exactly the momentum at which\n", + "[devos2022](@cite) locates the minimum of the dispersion, where that reference reports a gap\n", + "of $\\Delta / J = 0.0263$.\n", + "\n", + "The soft-mode *location* is reproduced cleanly, but the gap *value* is not, and it is worth\n", + "being precise about why.\n", + "The quasiparticle energies plotted above are a genuine variational upper bound on the\n", + "excitation energies only when they are measured relative to the *exact* ground state.\n", + "Here the ground state is a deliberately lightweight, modest-bond-dimension uniform $SU(3)$\n", + "MPS, and for this near-critical model it is not converged to that accuracy.\n", + "Because the reference state sits slightly too high in energy, the ansatz can find \"excitations\"\n", + "that lie below it, so near $k = 2\\pi/3$ the computed branch drops to essentially zero and even\n", + "slightly negative.\n", + "The dispersion here is therefore *not* a variational estimate of the true gap: the small value\n", + "$\\Delta / J = 0.0263$ is quoted from [devos2022](@cite) and would require a carefully\n", + "converged, large-bond-dimension calculation to reproduce.\n", + "What this example does show is the complete $SU(3)$-symmetric workflow — Hamiltonian\n", + "construction, VUMPS ground state, and a momentum-resolved excitation scan — with the soft mode\n", + "correctly located at $k = 2\\pi/3$." + ] + }, + { + "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.12.6", + "language": "julia", + "name": "julia-1.12" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.12.6" + } + }, + "nbformat": 4, + "nbformat_minor": 3 +} \ No newline at end of file diff --git a/docs/src/examples/groundstates/0.tfim-groundstate/figure-1.png b/docs/src/examples/groundstates/0.tfim-groundstate/figure-1.png new file mode 100644 index 000000000..bfa456242 Binary files /dev/null and b/docs/src/examples/groundstates/0.tfim-groundstate/figure-1.png differ diff --git a/docs/src/examples/groundstates/0.tfim-groundstate/figure-2.png b/docs/src/examples/groundstates/0.tfim-groundstate/figure-2.png new file mode 100644 index 000000000..19e71c3ea Binary files /dev/null and b/docs/src/examples/groundstates/0.tfim-groundstate/figure-2.png differ diff --git a/docs/src/examples/groundstates/0.tfim-groundstate/figure-3.png b/docs/src/examples/groundstates/0.tfim-groundstate/figure-3.png new file mode 100644 index 000000000..31b1c357a Binary files /dev/null and b/docs/src/examples/groundstates/0.tfim-groundstate/figure-3.png differ diff --git a/docs/src/examples/groundstates/0.tfim-groundstate/index.md b/docs/src/examples/groundstates/0.tfim-groundstate/index.md new file mode 100644 index 000000000..9552cde5e --- /dev/null +++ b/docs/src/examples/groundstates/0.tfim-groundstate/index.md @@ -0,0 +1,199 @@ +```@meta +EditURL = "../../../../../examples/groundstates/0.tfim-groundstate/main.jl" +``` + +[![](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/QuantumKitHub/MPSKit.jl/gh-pages?filepath=dev/examples/groundstates/0.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/0.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/0.tfim-groundstate) + +# The transverse-field Ising model: a complete ground-state study + +This example is the bridge from the introductory tutorials into the research-grade +gallery. +If you have worked through [Your first ground state](@ref tutorial_first_groundstate) and +[The thermodynamic limit](@ref tutorial_thermodynamic_limit) you already know every +individual tool used here; the goal now is to *assemble* them into one coherent case +study of a genuine quantum phase transition. + +We use the same transverse-field Ising model (TFIM) as the tutorials, on a chain of +spin-1/2 sites: + +```math +H = -J\left(\sum_{\langle i,j\rangle} \sigma^z_i \sigma^z_j + g\sum_i \sigma^x_i\right), +``` + +where the first sum runs over neighbouring pairs, ``J`` sets the energy scale, and the +dimensionless field ``g`` tunes the competition between the ``\sigma^z\sigma^z`` +interaction and the transverse ``\sigma^x`` field. +The model has a quantum critical point at ``g = 1``. + +Rather than looking at a single field value, we will scan ``g`` across the transition and +diagnose it three independent ways, comparing a *finite* chain against a calculation +performed *directly in the thermodynamic limit*: + +1. the order parameter ``|\langle\sigma^z\rangle|``, computed both for a finite chain and + for an infinite chain, in one figure; +2. the entanglement entropy of the infinite state; +3. the correlation length of the infinite state. + +All three should point at the same place — that agreement is the payoff. + +We take the model and lattice from MPSKitModels, the tensor backend from TensorKit, and +Plots for the figures. The Pauli operators `σᶻ`, `σˣ` are re-exported by MPSKitModels. + +````julia +using MPSKit, MPSKitModels, TensorKit, Plots +```` + +## Shared parameters + +We fix a finite chain length `L`, a bond dimension `D` (the accuracy knob, see +[Controlling bond dimension](@ref howto_bond_dimension)), and the set of field values to +scan. +`D` is kept modest so the whole page runs in a couple of minutes; increasing it sharpens +the infinite-state diagnostics below (the finite-chain curve responds to `D` in a less +obvious way, as we will see). + +````julia +L = 16 +D = 8 +g_values = 0.1:0.1:2.0 +```` + +```` +0.1:0.1:2.0 +```` + +## 1. Finite versus infinite magnetization + +We compute the order parameter ``|\langle\sigma^z\rangle|`` two ways at every field value. + +For the **finite** calculation we use an open chain of `L` sites, exactly as in the +tutorial, and optimize with [`DMRG`](@ref). +We average ``\langle\sigma^z_i\rangle`` over the sites and take the absolute value: the +exact finite-`L` ground state is symmetric, but DMRG lands on one of the two +symmetry-broken states with an arbitrary sign (see the discussion in +[Your first ground state](@ref tutorial_first_groundstate)). + +````julia +ψ₀_finite = FiniteMPS(L, ℂ^2, ℂ^D) +M_finite = map(g_values) do g + H = transverse_field_ising(FiniteChain(L); g = g) + ψ, = find_groundstate(ψ₀_finite, H, DMRG(; verbosity = 0)) + return abs(sum(expectation_value(ψ, i => σᶻ()) for i in 1:L)) / L +end; +```` + +For the **infinite** calculation we drop the lattice argument to build the Hamiltonian on +the infinite chain, use an [`InfiniteMPS`](@ref), and optimize with [`VUMPS`](@ref). +We keep every optimized infinite state, because we will reuse them for the entropy and +correlation-length diagnostics below. + +````julia +ψ₀_infinite = InfiniteMPS(ℂ^2, ℂ^D) +states_infinite = map(g_values) do g + H = transverse_field_ising(; g = g) + ψ, = find_groundstate(ψ₀_infinite, H, VUMPS(; verbosity = 0)) + return ψ +end; +```` + +The order parameter of a translation-invariant state is just ``\langle\sigma^z\rangle`` on +a single site of the unit cell; we again take the absolute value, because on the ordered +side the infinite state settles into one of the two symmetry-broken ground states (see +[The thermodynamic limit](@ref tutorial_thermodynamic_limit)). + +````julia +M_infinite = [abs(expectation_value(ψ, 1 => σᶻ())) for ψ in states_infinite]; +```` + +Plotting both curves in a single figure lets us compare them directly. + +````julia +p_magnetization = plot(; + xlabel = "g", ylabel = "|⟨σᶻ⟩|", title = "TFIM order parameter", legend = :bottomleft +) +scatter!(p_magnetization, g_values, M_finite; label = "finite chain, L = $L, D = $D") +scatter!(p_magnetization, g_values, M_infinite; label = "infinite, D = $D") +vline!(p_magnetization, [1.0]; color = "gray", linestyle = :dash, label = "g = 1") +p_magnetization +```` + +![](figure-1.png) + +Both calculations agree deep in either phase, but near the transition they tell very different stories. +The infinite curve stays on its ordered branch essentially up to `g = 1` and then collapses: it locates the critical point cleanly. +The finite-chain curve instead drops to zero far earlier — at this `L` and `D` the variational optimum on the open chain switches from the symmetry-broken branch to the exactly symmetric ground state, whose magnetization vanishes. +Where that switch happens is set by `L` and `D`, not by the physics; the same sweep at `D = 4` in [Your first ground state](@ref tutorial_first_groundstate) puts it elsewhere. +That is the real lesson of this panel: the finite-chain order parameter is dominated by which state the algorithm selects, while the calculation performed directly in the thermodynamic limit pins the transition at `g = 1`. + +## 2. Entanglement entropy across the transition + +Entanglement is a hallmark of criticality: it is bounded away from the critical point but +grows sharply as we approach it. +For an [`InfiniteMPS`](@ref), [`entropy`](@ref) returns the von Neumann entanglement entropy +per bond, one value for each site of the unit cell. +Our unit cell has a single site, so we take the one entry with `only`. + +````julia +S_infinite = [real(only(entropy(ψ))) for ψ in states_infinite] +p_entropy = scatter( + g_values, S_infinite; + xlabel = "g", ylabel = "entanglement entropy S", title = "TFIM entanglement entropy", + legend = false +) +vline!(p_entropy, [1.0]; color = "gray", linestyle = :dash) +p_entropy +```` + +![](figure-2.png) + +The entropy peaks near `g = 1`. +That peak is the entanglement signature of the phase transition: at criticality +correlations become long-ranged and the ground state is at its most entangled, whereas deep +in either phase the state is closer to a simple product and the entropy is small. + +## 3. Correlation length across the transition + +The [`correlation_length`](@ref) measures how far apart two spins can still influence each +other; it is extracted from the transfer-matrix spectrum of the uniform infinite state and +has no finite-chain analogue. +It grows toward criticality, so we plot it on a logarithmic vertical axis to make the +growth visible. + +````julia +ξ_infinite = [correlation_length(ψ) for ψ in states_infinite] +p_xi = scatter( + g_values, ξ_infinite; + xlabel = "g", ylabel = "correlation length ξ", yscale = :log10, + title = "TFIM correlation length", legend = false +) +vline!(p_xi, [1.0]; color = "gray", linestyle = :dash) +p_xi +```` + +![](figure-3.png) + +The correlation length peaks near `g = 1` as well. +At a genuine critical point it would diverge, but a finite bond dimension `D` can only +capture correlations out to a finite range, so what we measure is large-but-capped rather +than infinite — the peak grows and sharpens as `D` is increased. + +## What you now have + +Three independent diagnostics — the order parameter, the entanglement entropy, and the +correlation length — all locate the transition of the transverse-field Ising model near +`g = 1`, and the finite-versus-infinite comparison shows concretely why the thermodynamic +limit is the right place to measure it. + +From here the gallery goes further. +The Ising CFT example extracts the momentum-resolved excitation spectrum right at +criticality and matches it to the predictions of conformal field theory, turning the "there +is a critical point near `g = 1`" of this page into a quantitative fingerprint of *which* +critical theory it is. +Every curve on this page also sharpens if you rerun it at a larger bond dimension `D`. + +--- + +*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* + diff --git a/docs/src/examples/groundstates/0.tfim-groundstate/main.ipynb b/docs/src/examples/groundstates/0.tfim-groundstate/main.ipynb new file mode 100644 index 000000000..12748dde4 --- /dev/null +++ b/docs/src/examples/groundstates/0.tfim-groundstate/main.ipynb @@ -0,0 +1,310 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# The transverse-field Ising model: a complete ground-state study\n", + "\n", + "This example is the bridge from the introductory tutorials into the research-grade\n", + "gallery.\n", + "If you have worked through Your first ground state and\n", + "The thermodynamic limit you already know every\n", + "individual tool used here; the goal now is to *assemble* them into one coherent case\n", + "study of a genuine quantum phase transition.\n", + "\n", + "We use the same transverse-field Ising model (TFIM) as the tutorials, on a chain of\n", + "spin-1/2 sites:\n", + "\n", + "$$\n", + "H = -J\\left(\\sum_{\\langle i,j\\rangle} \\sigma^z_i \\sigma^z_j + g\\sum_i \\sigma^x_i\\right),\n", + "$$\n", + "\n", + "where the first sum runs over neighbouring pairs, $J$ sets the energy scale, and the\n", + "dimensionless field $g$ tunes the competition between the $\\sigma^z\\sigma^z$\n", + "interaction and the transverse $\\sigma^x$ field.\n", + "The model has a quantum critical point at $g = 1$.\n", + "\n", + "Rather than looking at a single field value, we will scan $g$ across the transition and\n", + "diagnose it three independent ways, comparing a *finite* chain against a calculation\n", + "performed *directly in the thermodynamic limit*:\n", + "\n", + "1. the order parameter $|\\langle\\sigma^z\\rangle|$, computed both for a finite chain and\n", + " for an infinite chain, in one figure;\n", + "2. the entanglement entropy of the infinite state;\n", + "3. the correlation length of the infinite state.\n", + "\n", + "All three should point at the same place — that agreement is the payoff." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We take the model and lattice from MPSKitModels, the tensor backend from TensorKit, and\n", + "Plots for the figures. The Pauli operators `σᶻ`, `σˣ` are re-exported by MPSKitModels." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "using MPSKit, MPSKitModels, TensorKit, Plots" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Shared parameters\n", + "\n", + "We fix a finite chain length `L`, a bond dimension `D` (the accuracy knob, see\n", + "Controlling bond dimension), and the set of field values to\n", + "scan.\n", + "`D` is kept modest so the whole page runs in a couple of minutes; increasing it sharpens\n", + "the infinite-state diagnostics below (the finite-chain curve responds to `D` in a less\n", + "obvious way, as we will see)." + ] + }, + { + "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": [ + "## 1. Finite versus infinite magnetization\n", + "\n", + "We compute the order parameter $|\\langle\\sigma^z\\rangle|$ two ways at every field value.\n", + "\n", + "For the **finite** calculation we use an open chain of `L` sites, exactly as in the\n", + "tutorial, and optimize with `DMRG`.\n", + "We average $\\langle\\sigma^z_i\\rangle$ over the sites and take the absolute value: the\n", + "exact finite-`L` ground state is symmetric, but DMRG lands on one of the two\n", + "symmetry-broken states with an arbitrary sign (see the discussion in\n", + "Your first ground state)." + ] + }, + { + "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(FiniteChain(L); g = g)\n", + " ψ, = find_groundstate(ψ₀_finite, H, DMRG(; verbosity = 0))\n", + " return abs(sum(expectation_value(ψ, i => σᶻ()) for i in 1:L)) / L\n", + "end;" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For the **infinite** calculation we drop the lattice argument to build the Hamiltonian on\n", + "the infinite chain, use an `InfiniteMPS`, and optimize with `VUMPS`.\n", + "We keep every optimized infinite state, because we will reuse them for the entropy and\n", + "correlation-length diagnostics below." + ] + }, + { + "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 = g)\n", + " ψ, = find_groundstate(ψ₀_infinite, H, VUMPS(; verbosity = 0))\n", + " return ψ\n", + "end;" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The order parameter of a translation-invariant state is just $\\langle\\sigma^z\\rangle$ on\n", + "a single site of the unit cell; we again take the absolute value, because on the ordered\n", + "side the infinite state settles into one of the two symmetry-broken ground states (see\n", + "The thermodynamic limit)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "M_infinite = [abs(expectation_value(ψ, 1 => σᶻ())) for ψ in states_infinite];" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plotting both curves in a single figure lets us compare them directly." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "p_magnetization = plot(;\n", + " xlabel = \"g\", ylabel = \"|⟨σᶻ⟩|\", title = \"TFIM order parameter\", legend = :bottomleft\n", + ")\n", + "scatter!(p_magnetization, g_values, M_finite; label = \"finite chain, L = $L, D = $D\")\n", + "scatter!(p_magnetization, g_values, M_infinite; label = \"infinite, D = $D\")\n", + "vline!(p_magnetization, [1.0]; color = \"gray\", linestyle = :dash, label = \"g = 1\")\n", + "p_magnetization" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Both calculations agree deep in either phase, but near the transition they tell very different stories.\n", + "The infinite curve stays on its ordered branch essentially up to `g = 1` and then collapses: it locates the critical point cleanly.\n", + "The finite-chain curve instead drops to zero far earlier — at this `L` and `D` the variational optimum on the open chain switches from the symmetry-broken branch to the exactly symmetric ground state, whose magnetization vanishes.\n", + "Where that switch happens is set by `L` and `D`, not by the physics; the same sweep at `D = 4` in Your first ground state puts it elsewhere.\n", + "That is the real lesson of this panel: the finite-chain order parameter is dominated by which state the algorithm selects, while the calculation performed directly in the thermodynamic limit pins the transition at `g = 1`." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. Entanglement entropy across the transition\n", + "\n", + "Entanglement is a hallmark of criticality: it is bounded away from the critical point but\n", + "grows sharply as we approach it.\n", + "For an `InfiniteMPS`, `entropy` returns the von Neumann entanglement entropy\n", + "per bond, one value for each site of the unit cell.\n", + "Our unit cell has a single site, so we take the one entry with `only`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "S_infinite = [real(only(entropy(ψ))) for ψ in states_infinite]\n", + "p_entropy = scatter(\n", + " g_values, S_infinite;\n", + " xlabel = \"g\", ylabel = \"entanglement entropy S\", title = \"TFIM entanglement entropy\",\n", + " legend = false\n", + ")\n", + "vline!(p_entropy, [1.0]; color = \"gray\", linestyle = :dash)\n", + "p_entropy" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The entropy peaks near `g = 1`.\n", + "That peak is the entanglement signature of the phase transition: at criticality\n", + "correlations become long-ranged and the ground state is at its most entangled, whereas deep\n", + "in either phase the state is closer to a simple product and the entropy is small." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 3. Correlation length across the transition\n", + "\n", + "The `correlation_length` measures how far apart two spins can still influence each\n", + "other; it is extracted from the transfer-matrix spectrum of the uniform infinite state and\n", + "has no finite-chain analogue.\n", + "It grows toward criticality, so we plot it on a logarithmic vertical axis to make the\n", + "growth visible." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ξ_infinite = [correlation_length(ψ) for ψ in states_infinite]\n", + "p_xi = scatter(\n", + " g_values, ξ_infinite;\n", + " xlabel = \"g\", ylabel = \"correlation length ξ\", yscale = :log10,\n", + " title = \"TFIM correlation length\", legend = false\n", + ")\n", + "vline!(p_xi, [1.0]; color = \"gray\", linestyle = :dash)\n", + "p_xi" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The correlation length peaks near `g = 1` as well.\n", + "At a genuine critical point it would diverge, but a finite bond dimension `D` can only\n", + "capture correlations out to a finite range, so what we measure is large-but-capped rather\n", + "than infinite — the peak grows and sharpens as `D` is increased." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## What you now have\n", + "\n", + "Three independent diagnostics — the order parameter, the entanglement entropy, and the\n", + "correlation length — all locate the transition of the transverse-field Ising model near\n", + "`g = 1`, and the finite-versus-infinite comparison shows concretely why the thermodynamic\n", + "limit is the right place to measure it.\n", + "\n", + "From here the gallery goes further.\n", + "The Ising CFT example extracts the momentum-resolved excitation spectrum right at\n", + "criticality and matches it to the predictions of conformal field theory, turning the \"there\n", + "is a critical point near `g = 1`\" of this page into a quantitative fingerprint of *which*\n", + "critical theory it is.\n", + "Every curve on this page also sharpens if you rerun it at a larger bond dimension `D`." + ] + }, + { + "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.12.6", + "language": "julia", + "name": "julia-1.12" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.12.6" + } + }, + "nbformat": 4, + "nbformat_minor": 3 +} \ No newline at end of file diff --git a/examples/Cache.toml b/examples/Cache.toml index f3e9fa9fe..d65c505c7 100644 --- a/examples/Cache.toml +++ b/examples/Cache.toml @@ -1,4 +1,5 @@ [excitations] +"1.su3-heisenberg" = "a57dbdf1b92f6c589c086b166dab64f4f224813ef632e1668c34c9607f005f9f" "0.haldane" = "c09df36c3be5cd452bec564d22f6e2815474a9dab71dd5528013a3fb46b2f7c0" [dynamics] @@ -13,4 +14,5 @@ "2.hubbard" = "29e2469ac9307f1242bdfb76b97c6a2e04dc087f5d4b6c45c8b7a15df12ad36c" "3.bose-hubbard" = "cf0d9a543e784dc6053e413780d19e1b59bf956dd4598752fdf664804dd68ce6" "1.xxz-heisenberg" = "dae5f29dcad5fcfaffd12a7847916adb92c955276ce24166fd3f9c9534c50784" +"0.tfim-groundstate" = "a87cc89caf8df47285e25f6cec4ff130ea6ed88d4f73d59f8e7dc2e799ece1b0" "5.ising-cft" = "ebcae1e501347cb00a46d0889df1355707be77cf962f9ee0d16416393492529b" diff --git a/examples/Project.toml b/examples/Project.toml index 2e29ddd9c..1d6ade97b 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -12,6 +12,7 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" 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" TensorOperations = "6aa20fa7-93e2-5fca-9bc0-fbd0db3c71a2" diff --git a/examples/excitations/1.su3-heisenberg/main.jl b/examples/excitations/1.su3-heisenberg/main.jl new file mode 100644 index 000000000..3c244f281 --- /dev/null +++ b/examples/excitations/1.su3-heisenberg/main.jl @@ -0,0 +1,190 @@ +md""" +# The SU(3) Heisenberg chain + +The spin-1 Heisenberg chain famously has a gapped, symmetry-protected topological ground +state: the Haldane phase. +A natural question is what happens when the ``SU(2)`` symmetry of the spins is promoted to a +larger ``SU(N)`` symmetry. +In this example we study the ``SU(3)`` generalization introduced in [devos2022](@cite): a +uniform chain where every site carries the fully-symmetric ``[3\,0\,0]`` irrep of ``SU(3)`` +(the ten-dimensional representation), coupled by a nearest-neighbour ``SU(3)``-invariant +Heisenberg interaction + +```math +H = J \sum_{\langle i, j \rangle} \sum_{a=1}^{8} T^a_i T^a_j , +``` + +where the ``T^a`` are the eight generators of ``SU(3)``. +This is the direct analog of the spin-1 chain, with the three generators of ``SU(2)`` +replaced by the eight generators of ``SU(3)``. +The goal of this example is to exercise the full ``SU(N)`` toolchain end to end — building the +symmetric Hamiltonian, optimizing a uniform ground state with VUMPS, and scanning an +excitation dispersion — rather than to pin down the model's (very small) gap to high accuracy; +we return to that point in the closing discussion. + +Working with the full non-abelian ``SU(3)`` symmetry is what makes this tractable. +Just as ``SU(2)`` symmetry drastically shrinks the spin-1 problem (see +[Symmetries](@ref concept_symmetries)), imposing ``SU(3)`` block-diagonalizes every tensor +along its irreducible representations, so that a modest number of variational parameters +already captures a large physical bond dimension. +The Clebsch–Gordan coefficients that this requires are provided by +[SUNRepresentations.jl](https://github.com/QuantumKitHub/SUNRepresentations.jl), whose +`SU3Irrep` sectors plug directly into TensorKit's graded vector spaces. + +!!! warning "First-run cost" + SUNRepresentations.jl computes ``SU(N)`` Clebsch–Gordan coefficients on demand and + caches them both in memory and on disk (via `Scratch.jl`). + The **first** time a given coefficient is needed the computation can be slow; later runs + reuse the on-disk cache. + On a machine without a persisted scratch cache — such as a fresh continuous-integration + runner — this example therefore pays the full coefficient-generation cost once. +""" + +using TensorKit +using SUNRepresentations +using MPSKit +using Plots + +# 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 obtain from its highest weight and +represent as a one-dimensional graded space. +""" + +I300 = SU3Irrep(3, 0, 0) +Vphys = Vect[SU3Irrep](I300 => 1) +@show dim(I300) + +md""" +To build the interaction we use the same Casimir trick as for the ``SU(2)`` chain. +Writing ``C_2 = \sum_a T^a T^a`` for the quadratic Casimir operator, the two-site coupling is +diagonal in the total ``SU(3)`` charge ``c`` of a bond: + +```math +\sum_a T^a_i T^a_j = \tfrac{1}{2}\left( C_2(c) - C_2([3\,0\,0]) - C_2([3\,0\,0]) \right). +``` + +SUNRepresentations.jl provides the quadratic Casimir through `casimir(2, irrep)`. +Its normalization is the standard one with ``\mathrm{tr}(T^a T^b) = \tfrac{1}{2}\delta^{ab}``: +the fundamental ``[1\,0\,0]`` has ``C_2 = 4/3 = (N^2-1)/(2N)``. +""" + +@show casimir(2, SU3Irrep(1, 0, 0)) +@show casimir(2, I300) + +md""" +The nearest-neighbour term is then a single TensorKit operator that acts as a scalar on each +fusion channel of ``[3\,0\,0] \otimes [3\,0\,0]``. +Constructing it by iterating over the blocks — one block per total charge ``c`` — mirrors the +``SU(2)`` construction, but now with `casimir(2, ⋅)` supplying the eigenvalues. +""" + +function su3_heisenberg(Vphys; J = 1.0) + site = SU3Irrep(3, 0, 0) + c2_site = casimir(2, site) + 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""" +We can read off the physics directly from the fusion channels. +The product ``[3\,0\,0] \otimes [3\,0\,0]`` decomposes into four irreps, and the +antiferromagnetic coupling (``J > 0``) favours the channels of lowest Casimir: +""" + +for (c, n) in directproduct(I300, I300) + coupling = (casimir(2, c) - 2 * casimir(2, I300)) / 2 + println(rpad(string(c), 22), " dim = ", rpad(dim(c), 4), " coupling = ", coupling) +end + +md""" +## Ground state + +The ``[3\,0\,0]`` irrep has trivial triality (congruency class ``0``), so — like the +integer-spin Haldane chain — a translation-invariant state on a single-site unit cell is +admissible. +Because the physical charge has triality ``0`` it cannot shift the triality of the virtual +bonds, so every virtual space must sit in a single congruency class; here we use the trivial +class ``0``, which we found to give the lowest variational energy. +Following [devos2022](@cite), we seed the optimization with a random ``SU(3)``-symmetric +uniform state and optimize with VUMPS. +The bond dimension is kept deliberately modest here so that the page builds quickly. +""" + +Vvirt = Vect[SU3Irrep]( + SU3Irrep(0, 0, 0) => 8, + SU3Irrep(2, 1, 0) => 6, + SU3Irrep(3, 0, 0) => 3, + SU3Irrep(3, 3, 0) => 3, + SU3Irrep(4, 2, 0) => 2, +) +@show dim(Vvirt) + +ψ₀ = InfiniteMPS([Vphys], [Vvirt]) +ψ, envs, δ = find_groundstate(ψ₀, H, VUMPS(; maxiter = 200, tol = 1.0e-10, verbosity = 1)) +E₀ = real(expectation_value(ψ, H)) +println("ground-state energy per site: E₀ = $E₀") + +md""" +## The excitation spectrum + +On top of the ground state we compute excitations with the quasiparticle ansatz. +In the thermodynamic limit each excitation carries a definite momentum ``k`` and a definite +``SU(3)`` charge. +Following [devos2022](@cite), we look in the adjoint ``[2\,1\,0]`` sector, where that reference +finds the lowest branch, and scan across half of the Brillouin zone (the other half follows by +reflection). +""" + +sector = SU3Irrep(2, 1, 0) +kspace = range(0, π, 10) +Es, _ = excitations(H, QuasiparticleAnsatz(), kspace, ψ, envs; sector, verbosity = 0) + +Δ, idx = findmin(real.(Es)) +println("branch minimum ΔE/J = $Δ at k = $(kspace[idx])") + +#+ + +plot( + kspace, real.(Es); + xaxis = "momentum k", yaxis = "ΔE / J", label = "[2 1 0] branch", + title = "SU(3) [3 0 0] excitation dispersion" +) +hline!([0.0]; color = :gray, linestyle = :dash, label = "") +vline!([2π / 3]; color = :gray, linestyle = :dot, label = "k = 2π/3") + +md""" +The branch has a pronounced soft mode at ``k = 2\pi/3``, exactly the momentum at which +[devos2022](@cite) locates the minimum of the dispersion, where that reference reports a gap +of ``\Delta / J = 0.0263``. + +The soft-mode *location* is reproduced cleanly, but the gap *value* is not, and it is worth +being precise about why. +The quasiparticle energies plotted above are a genuine variational upper bound on the +excitation energies only when they are measured relative to the *exact* ground state. +Here the ground state is a deliberately lightweight, modest-bond-dimension uniform ``SU(3)`` +MPS, and for this near-critical model it is not converged to that accuracy. +Because the reference state sits slightly too high in energy, the ansatz can find "excitations" +that lie below it, so near ``k = 2\pi/3`` the computed branch drops to essentially zero and even +slightly negative. +The dispersion here is therefore *not* a variational estimate of the true gap: the small value +``\Delta / J = 0.0263`` is quoted from [devos2022](@cite) and would require a carefully +converged, large-bond-dimension calculation to reproduce. +What this example does show is the complete ``SU(3)``-symmetric workflow — Hamiltonian +construction, VUMPS ground state, and a momentum-resolved excitation scan — with the soft mode +correctly located at ``k = 2\pi/3``. +""" diff --git a/examples/groundstates/0.tfim-groundstate/main.jl b/examples/groundstates/0.tfim-groundstate/main.jl new file mode 100644 index 000000000..09b3ae7e5 --- /dev/null +++ b/examples/groundstates/0.tfim-groundstate/main.jl @@ -0,0 +1,185 @@ +md""" +# The transverse-field Ising model: a complete ground-state study + +This example is the bridge from the introductory tutorials into the research-grade +gallery. +If you have worked through [Your first ground state](@ref tutorial_first_groundstate) and +[The thermodynamic limit](@ref tutorial_thermodynamic_limit) you already know every +individual tool used here; the goal now is to *assemble* them into one coherent case +study of a genuine quantum phase transition. + +We use the same transverse-field Ising model (TFIM) as the tutorials, on a chain of +spin-1/2 sites: + +```math +H = -J\left(\sum_{\langle i,j\rangle} \sigma^z_i \sigma^z_j + g\sum_i \sigma^x_i\right), +``` + +where the first sum runs over neighbouring pairs, ``J`` sets the energy scale, and the +dimensionless field ``g`` tunes the competition between the ``\sigma^z\sigma^z`` +interaction and the transverse ``\sigma^x`` field. +The model has a quantum critical point at ``g = 1``. + +Rather than looking at a single field value, we will scan ``g`` across the transition and +diagnose it three independent ways, comparing a *finite* chain against a calculation +performed *directly in the thermodynamic limit*: + +1. the order parameter ``|\langle\sigma^z\rangle|``, computed both for a finite chain and + for an infinite chain, in one figure; +2. the entanglement entropy of the infinite state; +3. the correlation length of the infinite state. + +All three should point at the same place — that agreement is the payoff. +""" + +# We take the model and lattice from MPSKitModels, the tensor backend from TensorKit, and +# Plots for the figures. The Pauli operators `σᶻ`, `σˣ` are re-exported by MPSKitModels. + +using MPSKit, MPSKitModels, TensorKit, Plots + +md""" +## Shared parameters + +We fix a finite chain length `L`, a bond dimension `D` (the accuracy knob, see +[Controlling bond dimension](@ref howto_bond_dimension)), and the set of field values to +scan. +`D` is kept modest so the whole page runs in a couple of minutes; increasing it sharpens +the infinite-state diagnostics below (the finite-chain curve responds to `D` in a less +obvious way, as we will see). +""" + +L = 16 +D = 8 +g_values = 0.1:0.1:2.0 + +md""" +## 1. Finite versus infinite magnetization + +We compute the order parameter ``|\langle\sigma^z\rangle|`` two ways at every field value. + +For the **finite** calculation we use an open chain of `L` sites, exactly as in the +tutorial, and optimize with [`DMRG`](@ref). +We average ``\langle\sigma^z_i\rangle`` over the sites and take the absolute value: the +exact finite-`L` ground state is symmetric, but DMRG lands on one of the two +symmetry-broken states with an arbitrary sign (see the discussion in +[Your first ground state](@ref tutorial_first_groundstate)). +""" + +ψ₀_finite = FiniteMPS(L, ℂ^2, ℂ^D) +M_finite = map(g_values) do g + H = transverse_field_ising(FiniteChain(L); g = g) + ψ, = find_groundstate(ψ₀_finite, H, DMRG(; verbosity = 0)) + return abs(sum(expectation_value(ψ, i => σᶻ()) for i in 1:L)) / L +end; + +md""" +For the **infinite** calculation we drop the lattice argument to build the Hamiltonian on +the infinite chain, use an [`InfiniteMPS`](@ref), and optimize with [`VUMPS`](@ref). +We keep every optimized infinite state, because we will reuse them for the entropy and +correlation-length diagnostics below. +""" + +ψ₀_infinite = InfiniteMPS(ℂ^2, ℂ^D) +states_infinite = map(g_values) do g + H = transverse_field_ising(; g = g) + ψ, = find_groundstate(ψ₀_infinite, H, VUMPS(; verbosity = 0)) + return ψ +end; + +md""" +The order parameter of a translation-invariant state is just ``\langle\sigma^z\rangle`` on +a single site of the unit cell; we again take the absolute value, because on the ordered +side the infinite state settles into one of the two symmetry-broken ground states (see +[The thermodynamic limit](@ref tutorial_thermodynamic_limit)). +""" + +M_infinite = [abs(expectation_value(ψ, 1 => σᶻ())) for ψ in states_infinite]; + +md""" +Plotting both curves in a single figure lets us compare them directly. +""" + +p_magnetization = plot(; + xlabel = "g", ylabel = "|⟨σᶻ⟩|", title = "TFIM order parameter", legend = :bottomleft +) +scatter!(p_magnetization, g_values, M_finite; label = "finite chain, L = $L, D = $D") +scatter!(p_magnetization, g_values, M_infinite; label = "infinite, D = $D") +vline!(p_magnetization, [1.0]; color = "gray", linestyle = :dash, label = "g = 1") +p_magnetization + +md""" +Both calculations agree deep in either phase, but near the transition they tell very different stories. +The infinite curve stays on its ordered branch essentially up to `g = 1` and then collapses: it locates the critical point cleanly. +The finite-chain curve instead drops to zero far earlier — at this `L` and `D` the variational optimum on the open chain switches from the symmetry-broken branch to the exactly symmetric ground state, whose magnetization vanishes. +Where that switch happens is set by `L` and `D`, not by the physics; the same sweep at `D = 4` in [Your first ground state](@ref tutorial_first_groundstate) puts it elsewhere. +That is the real lesson of this panel: the finite-chain order parameter is dominated by which state the algorithm selects, while the calculation performed directly in the thermodynamic limit pins the transition at `g = 1`. + +""" + +md""" +## 2. Entanglement entropy across the transition + +Entanglement is a hallmark of criticality: it is bounded away from the critical point but +grows sharply as we approach it. +For an [`InfiniteMPS`](@ref), [`entropy`](@ref) returns the von Neumann entanglement entropy +per bond, one value for each site of the unit cell. +Our unit cell has a single site, so we take the one entry with `only`. +""" + +S_infinite = [real(only(entropy(ψ))) for ψ in states_infinite] +p_entropy = scatter( + g_values, S_infinite; + xlabel = "g", ylabel = "entanglement entropy S", title = "TFIM entanglement entropy", + legend = false +) +vline!(p_entropy, [1.0]; color = "gray", linestyle = :dash) +p_entropy + +md""" +The entropy peaks near `g = 1`. +That peak is the entanglement signature of the phase transition: at criticality +correlations become long-ranged and the ground state is at its most entangled, whereas deep +in either phase the state is closer to a simple product and the entropy is small. +""" + +md""" +## 3. Correlation length across the transition + +The [`correlation_length`](@ref) measures how far apart two spins can still influence each +other; it is extracted from the transfer-matrix spectrum of the uniform infinite state and +has no finite-chain analogue. +It grows toward criticality, so we plot it on a logarithmic vertical axis to make the +growth visible. +""" + +ξ_infinite = [correlation_length(ψ) for ψ in states_infinite] +p_xi = scatter( + g_values, ξ_infinite; + xlabel = "g", ylabel = "correlation length ξ", yscale = :log10, + title = "TFIM correlation length", legend = false +) +vline!(p_xi, [1.0]; color = "gray", linestyle = :dash) +p_xi + +md""" +The correlation length peaks near `g = 1` as well. +At a genuine critical point it would diverge, but a finite bond dimension `D` can only +capture correlations out to a finite range, so what we measure is large-but-capped rather +than infinite — the peak grows and sharpens as `D` is increased. +""" + +md""" +## What you now have + +Three independent diagnostics — the order parameter, the entanglement entropy, and the +correlation length — all locate the transition of the transverse-field Ising model near +`g = 1`, and the finite-versus-infinite comparison shows concretely why the thermodynamic +limit is the right place to measure it. + +From here the gallery goes further. +The Ising CFT example extracts the momentum-resolved excitation spectrum right at +criticality and matches it to the predictions of conformal field theory, turning the "there +is a critical point near `g = 1`" of this page into a quantitative fingerprint of *which* +critical theory it is. +Every curve on this page also sharpens if you rerun it at a larger bond dimension `D`. +"""