From d6c8b2101fce0cf606285d4a78a64af3e9917513 Mon Sep 17 00:00:00 2001 From: Akshay Shankar Date: Fri, 31 Oct 2025 15:13:34 +0100 Subject: [PATCH 1/8] add bose hubbard example --- .../quantum1d/8.bose-hubbard/index.md | 346 ++++++++++++++ .../quantum1d/8.bose-hubbard/main.ipynb | 424 ++++++++++++++++++ examples/Cache.toml | 1 + examples/Project.toml | 1 + examples/quantum1d/8.bose-hubbard/main.jl | 265 +++++++++++ 5 files changed, 1037 insertions(+) create mode 100644 docs/src/examples/quantum1d/8.bose-hubbard/index.md create mode 100644 docs/src/examples/quantum1d/8.bose-hubbard/main.ipynb create mode 100644 examples/quantum1d/8.bose-hubbard/main.jl diff --git a/docs/src/examples/quantum1d/8.bose-hubbard/index.md b/docs/src/examples/quantum1d/8.bose-hubbard/index.md new file mode 100644 index 000000000..1e1a1af2e --- /dev/null +++ b/docs/src/examples/quantum1d/8.bose-hubbard/index.md @@ -0,0 +1,346 @@ +```@meta +EditURL = "../../../../../examples/quantum1d/8.bose-hubbard/main.jl" +``` + +[![](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/QuantumKitHub/MPSKit.jl/gh-pages?filepath=dev/examples/quantum1d/8.bose-hubbard/main.ipynb) +[![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](https://nbviewer.jupyter.org/github/QuantumKitHub/MPSKit.jl/blob/gh-pages/dev/examples/quantum1d/8.bose-hubbard/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/quantum1d/8.bose-hubbard) + +````julia +using Markdown +using MPSKit, MPSKitModels, TensorKit +using Plots, LaTeXStrings + +theme(:wong) +default(fontfamily="Computer Modern", label=nothing, dpi=100, framestyle=:box) +```` + +# 1D Bose-Hubbard model + +In this tutorial, we will explore the physics of the one-dimensional Bose–Hubbard model using matrix product states. For the most part, we replicate the results presented in [**Phys. Rev. B 105, 134502**](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.105.134502), which can be consulted for any statements in this tutorial that are not otherwise cited. The Hamiltonian under study is now defined as follows: + +$$H = -t \sum_{i} (\hat{a}_i^{\dagger} \hat{a}_{i+1} + \hat{a}_{i+1}^{\dagger} \hat{a}_i) + \frac{U}{2} \sum_i \hat{n}_i(\hat{n}_i - 1) - \mu \sum_i \hat{n}_i$$ + +where the bosonic creation and annihilation operators satisfy the canonical commutation relations (CCR): + +$$[\hat{a}_i, \hat{a}_j^{\dagger}] = \delta_{ij}.$$ + +Each lattice site hosts a local Hilbert space corresponding to bosonic occupation states $|n\rangle$, where $(n = 0, 1, 2, \ldots)$. Since this space is formally infinite-dimensional, numerical simulations typically impose a truncation at some maximum occupation number $(n_{\text{max}})$. Such a treatment is justified since it can be observed that the simulation results quickly converge with the cutoff if the filling fraction is kept sufficiently low. + +Within this truncated space, the local creation and annihilation operators are represented by finite-dimensional matrices. For example, with cutoff $n_{\text{max}}$, the annihilation operator takes the form + +$$ +\hat{a} = +\begin{bmatrix} +0 & \sqrt{1} & 0 & 0 & \cdots & 0 \\ +0 & 0 & \sqrt{2} & 0 & \cdots & 0 \\ +0 & 0 & 0 & \sqrt{3} & \cdots & 0 \\ +\vdots & & & \ddots & \ddots & \vdots \\ +0 & 0 & 0 & \cdots & 0 & \sqrt{n_{\text{max}}} \\ +0 & 0 & 0 & \cdots & 0 & 0 +\end{bmatrix} +$$ + +and the creation operator is simply its Hermitian conjugate, + +$$ +\hat{a}^\dagger = +\begin{bmatrix} +0 & 0 & 0 & \cdots & 0 & 0 \\ +\sqrt{1} & 0 & 0 & \cdots & 0 & 0 \\ +0 & \sqrt{2} & 0 & \cdots & 0 & 0 \\ +\vdots & & \ddots & \ddots & & \vdots \\ +0 & 0 & 0 & \cdots & 0 & 0 \\ +0 & 0 & 0 & \cdots & \sqrt{n_{\text{max}}} & 0 +\end{bmatrix} +$$ + +The number operator is then given by + +$$\hat{n} = \hat{a}^\dagger \hat{a} = \mathrm{diag}(0, 1, 2, \ldots, n_{\text{max}}).$$ + +Before moving on, notice that the Hamiltonian is uniform and translationally invariant. Typically, such models are studied on a finite chain of $N$ sites with periodic boundary conditions, but this introduces finite-size effects that are rather annoying to deal with. In contrast, the MPS framework allows us to work directly in the thermodynamic limit, avoiding such artifacts. We will follow this line of exploration in this tutorial and leave finite systems for another time. + +In order to work in the thermodynamic limit, we will have to create an `InfiniteMPS`. A complete specification of the MPS requires us to define the physical space and the virtual space of the constituent tensors. At this point is it useful to note that `MPSKit.jl` is powered by `TensorKit.jl` under the hood and has some very generic interfaces in order to allow imposing symmetries of all kinds. As a result, it is sometimes necessary to be a bit more explicit about what we want to do in terms of the vector spaces involved. In this case, we will not consider any symmetries and simply take the most naive approach of working within the `Trivial` sector. The physical space is then `ℂ^(nmax+1)`, and the virtual space is `ℂ^D` where $D$ is some integer chosen to be the bond dimension of the MPS, and `ℂ` is an alias for `ComplexSpace` (typeset as `\bbC`). As $D$ is increased, one increases the amount of entanglement, i.e, quantum correlations that can be captured by the state. + +````julia +cutoff, D = 4, 5 +initial_state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) +```` + +```` +single site InfiniteMPS: +│ ⋮ +│ C[1]: TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 1, Vector{ComplexF64}}(ComplexF64[0.97751+0.0im, 0.174631+0.0125955im, 0.104701+0.000310339im, 0.0116397-0.00627461im, 0.0464636+0.000250548im, 0.0+0.0im, 0.0143524+0.0im, 0.00360949+0.00128115im, -0.00368336-0.0047228im, 0.000692314+0.00211969im, 0.0+0.0im, 0.0+0.0im, 0.00655329+0.0im, 0.00552591+0.00129082im, -0.00303456+0.00143664im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0102319+0.0im, -0.00205334+0.00505879im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.00640026+0.0im], ℂ^5 ← ℂ^5) +├── AL[1]: TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[0.390809+0.214646im, 0.0489818+0.0751624im, 0.0497281+0.0204177im, 0.0263021-0.000245521im, 0.0120706+0.0101434im, 0.211349+0.301151im, 0.0250289+0.0560838im, 0.0305354+0.027863im, 0.0167841+3.42096e-5im, 0.0149195+0.010686im, 0.324269+0.45im, 0.0444069+0.0705503im, 0.0396325+0.0521653im, 0.00128175+0.00922583im, 0.0179868+0.0194204im, 0.29577+0.274281im, 0.0536787+0.0464249im, 0.0358399+0.0188768im, 0.0105704-0.00635067im, 0.00311747+0.0202148im, 0.234359+0.307432im, 0.0459379+0.0535219im, 0.0237564+0.0496054im, -0.00216268+0.0275897im, 0.00395079+0.000283255im, -0.24849+0.340759im, -0.0387098-0.00943688im, -0.066458+0.00845679im, -0.0827909-0.00877386im, 0.0442285+0.0196632im, 0.455119-0.222154im, 0.187925-0.0124478im, 0.0290096-0.0172291im, -0.0627177-0.0544896im, 0.0249878+0.0111796im, -0.258925-0.366019im, -0.0495616+0.00700246im, -0.0334674-0.0427697im, -0.0235555+0.00522338im, 0.0129344-0.023551im, 0.370142+0.3491im, 0.0503108+0.0969619im, 0.0315893+0.0386184im, 0.00266811+0.00258244im, 0.0500739+0.0254822im, 0.0565483+0.077467im, 0.0203963+0.0620186im, -0.00282502-0.0492105im, 0.0218008-0.057548im, -0.00380262+0.0459115im, 0.216276+0.596959im, 0.103859+0.0332121im, 0.031629+0.105885im, 0.0856072+0.0532032im, -0.0631615+0.0551744im, 0.0300362-0.137117im, -0.086399-0.0116169im, 0.0315047+0.0466518im, 0.01568+0.0662402im, 0.000724708-0.0406109im, 0.234667-0.308019im, 0.167079-0.192054im, 0.00466488-0.0928214im, 0.00539646-0.0186621im, 0.00962596-0.0344942im, -0.456697-0.139859im, -0.0242026-0.0025607im, -0.038756-0.014274im, 0.0144336+0.0956192im, -0.0437668-0.0752814im, 0.107807+0.107598im, -0.00732072+0.0396839im, 0.102896+0.0301685im, 0.0688301-0.0658399im, 0.0317364-0.00169938im, 0.0520469-0.259208im, -0.0400713-0.0759im, -0.0490438+0.0512775im, -0.141529+0.0757342im, 0.0237932-0.0950623im, -0.109863-0.118982im, 0.0226764-0.193153im, -0.094942-0.070543im, -0.0890496-0.0436445im, 0.0274223+0.0547519im, -0.039393-0.218089im, 0.0163183+0.21468im, 0.0132428-0.0222767im, 0.0425839-0.114359im, 0.0249832+0.0736523im, -0.121223-0.0674908im, -0.0640196+0.0505719im, -0.0662815+0.026577im, -0.0533572-0.126597im, 0.0568246+0.0547217im, 0.169584+0.743324im, -0.00516459+0.0778608im, -0.0791716+0.0630091im, -0.0966301+0.00510259im, 0.0637203+0.0805271im, 0.140845+0.171094im, 0.021739-0.177908im, -0.0762735+0.11312im, -0.191275+0.0465106im, 0.0304639-0.0648004im, 0.229036+0.518716im, 0.0905234+0.039573im, -0.0974247-0.0245568im, 0.0249835-0.092148im, -0.0546524+0.101555im, -0.222658-0.0990489im, -0.198631+0.169606im, -0.0657467+0.0130075im, 0.0849531-0.0726054im, -0.0655968+0.0338763im, -0.0709195-0.353206im, -0.119811-0.0898848im, -0.0932307+0.106653im, -0.10785-0.0538011im, 0.0929571+0.0313938im, 0.132025-0.162429im, 0.0960656-0.142876im, -0.0589741-0.111375im, -0.14054-0.142388im, 0.0681092+0.109553im], (ℂ^5 ⊗ ℂ^5) ← ℂ^5) +│ ⋮ + +```` + +This simply initializes a tensor filled with random entries (check out the documentation for other useful constructors). Next, we need the creation and annihilation operators. While we could construct them from scratch, here we will use `MPSKitModels.jl` instead that has predefined operators and models for most well-known lattice models. + +````julia +a_op = a_min(cutoff=cutoff) # creates a bosonic annihilation operator without any symmetries +display(a_op[]) +display((a_op'*a_op)[]) +```` + +The [] accessor lets us see the underlying array and indeed the operators are exactly what we require. Similarly, the Bose Hubbard model is also predefined (although we will construct our own variant later on). + +````julia +hamiltonian = bose_hubbard_model(InfiniteChain(1), cutoff=cutoff, U=1, mu=0.5, t=0.2) # It is not strictly required to pass InfiniteChain() and is only included for clarity; one may instead pass FiniteChain(N) as well +```` + +```` +single site InfiniteMPOHamiltonian{MPSKit.JordanMPOTensor{ComplexF64, TensorKit.ComplexSpace, Union{TensorKit.BraidingTensor{ComplexF64, TensorKit.ComplexSpace}, TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 2, Vector{ComplexF64}}}, TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}, TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 2, Vector{ComplexF64}}, TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 1, Vector{ComplexF64}}}}: +╷ ⋮ +┼ W[1]: JordanMPOTensor{ComplexF64, ComplexSpace, Union{BraidingTensor{ComplexF64, ComplexSpace}, TensorMap{ComplexF64, ComplexSpace, 2, 2, Vector{ComplexF64}}}, TensorMap{ComplexF64, ComplexSpace, 2, 1, Vector{ComplexF64}}, TensorMap{ComplexF64, ComplexSpace, 1, 2, Vector{ComplexF64}}, TensorMap{ComplexF64, ComplexSpace, 1, 1, Vector{ComplexF64}}}(((ℂ^1 ⊞ ℂ^2 ⊞ ℂ^1) ⊗ ⊞(ℂ^5)) ← (⊞(ℂ^5) ⊗ (ℂ^1 ⊞ ℂ^2 ⊞ ℂ^1)), SparseBlockTensorMap{Union{BraidingTensor{ComplexF64, ComplexSpace}, TensorMap{ComplexF64, ComplexSpace, 2, 2, Vector{ComplexF64}}}, ComplexF64, ComplexSpace, 2, 2, 4}(Dict{CartesianIndex{4}, Union{BraidingTensor{ComplexF64, ComplexSpace}, TensorMap{ComplexF64, ComplexSpace, 2, 2, Vector{ComplexF64}}}}(), (⊞(ℂ^2) ⊗ ⊞(ℂ^5)) ← (⊞(ℂ^5) ⊗ ⊞(ℂ^2))), SparseBlockTensorMap{TensorMap{ComplexF64, ComplexSpace, 2, 1, Vector{ComplexF64}}, ComplexF64, ComplexSpace, 2, 1, 3}(Dict{CartesianIndex{3}, TensorMap{ComplexF64, ComplexSpace, 2, 1, Vector{ComplexF64}}}(CartesianIndex(1, 1, 1)=>TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 3.16228+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 3.16228+0.0im, 9.66041e-16+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 4.47214+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, -5.54668e-31+0.0im, 4.47214+0.0im, -2.08077e-16+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 5.47723+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 5.47723+0.0im, -7.40208e-16+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 6.32456+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 6.32456+0.0im, 3.05151e-16+0.0im, 0.0+0.0im, 0.0+0.0im], (ℂ^2 ⊗ ℂ^5) ← ℂ^5)), (⊞(ℂ^2) ⊗ ⊞(ℂ^5)) ← ⊞(ℂ^5)), SparseBlockTensorMap{TensorMap{ComplexF64, ComplexSpace, 1, 2, Vector{ComplexF64}}, ComplexF64, ComplexSpace, 1, 2, 3}(Dict{CartesianIndex{3}, TensorMap{ComplexF64, ComplexSpace, 1, 2, Vector{ComplexF64}}}(CartesianIndex(1, 1, 1)=>TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 2, Vector{ComplexF64}}(ComplexF64[0.0+0.0im, -0.0632456+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, -0.0894427+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, -0.109545+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, -0.126491+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, -3.40958e-18+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, -0.0632456+0.0im, 0.0+0.0im, -3.19295e-18+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, -0.0894427+0.0im, 0.0+0.0im, -1.11022e-17+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, -0.109545+0.0im, 0.0+0.0im, 5.55112e-18+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, -0.126491+0.0im, 0.0+0.0im], ℂ^5 ← (ℂ^5 ⊗ ℂ^2))), ⊞(ℂ^5) ← (⊞(ℂ^5) ⊗ ⊞(ℂ^2))), SparseBlockTensorMap{TensorMap{ComplexF64, ComplexSpace, 1, 1, Vector{ComplexF64}}, ComplexF64, ComplexSpace, 1, 1, 2}(Dict{CartesianIndex{2}, TensorMap{ComplexF64, ComplexSpace, 1, 1, Vector{ComplexF64}}}(CartesianIndex(1, 1)=>TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 1, Vector{ComplexF64}}(ComplexF64[0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, -0.5+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 1.5+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 4.0+0.0im], ℂ^5 ← ℂ^5)), ⊞(ℂ^5) ← ⊞(ℂ^5))) +╵ ⋮ + +```` + +This has created the Hamiltonian operator as a matrix product operator (MPO) which is a convenient form to use in conjunction with MPS. Finally, the ground state optimization may be performed with either `iDMRG` or `VUMPS`. Both should take similar arguments but it is known that VUMPS is typically more efficient for these systems so we proceed with that. + +````julia +ground_state, _, _ = find_groundstate(initial_state, hamiltonian, VUMPS(tol=1e-6, verbosity=2, maxiter=200)) +println("Energy: ", expectation_value(ground_state, hamiltonian)) +```` + +```` +[ Info: VUMPS init: obj = +3.785371201331e-01 err = 5.9503e-01 +[ Info: VUMPS conv 48: obj = -6.757777651160e-01 err = 8.4647387076e-07 time = 13.78 sec +Energy: -0.6757777651160398 - 5.389840049968564e-17im + +```` + +This automatically runs the algorithm until a certain error [measure](https://github.com/QuantumKitHub/MPSKit.jl/blob/main/src/algorithms/toolbox.jl#L42) falls below the specified tolerance or the maximum iterations is reached. Let us wrap all this into a convenient function. + +````julia +function get_ground_state(mu, t, cutoff, D; kwargs...) + hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff=cutoff, U=1, mu=mu, t=t) + state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) + state, _, _ = find_groundstate(state, hamiltonian, VUMPS(; kwargs...)) + + return state +end + +ground_state = get_ground_state(0.5, 0.01, cutoff, D; tol=1e-6, verbosity=2, maxiter=500) +```` + +```` +single site InfiniteMPS: +│ ⋮ +│ C[1]: TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 1, Vector{ComplexF64}}(ComplexF64[0.54949+0.0im, 0.369229+0.0647129im, 0.178404+0.141066im, -0.315274+0.540671im, -0.123933+0.313108im, 0.0+0.0im, 0.0146972+0.0im, 0.00310867+0.00509383im, -0.00170132+0.0120718im, -0.00584434+0.00729265im, 0.0+0.0im, 0.0+0.0im, 0.0120375+0.0im, -0.00309885-0.000910003im, -0.00472919-0.00552589im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 3.69614e-6+0.0im, 1.57727e-6+3.97113e-7im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 3.35037e-6+0.0im], ℂ^5 ← ℂ^5) +├── AL[1]: TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[0.334572+0.0477944im, 0.0555237+0.176413im, 0.166278+0.112267im, -0.33895+0.132832im, -0.124122-0.025603im, 0.243524+0.200797im, 0.12627+0.147902im, 0.0280863+0.120709im, -0.321951+0.11476im, -0.154783+0.0931863im, 0.188335+0.150819im, 0.0985741+0.105361im, 0.0296448+0.0819545im, -0.241768+0.094309im, -0.120223+0.0771641im, 0.124161+0.204734im, 0.0575478+0.149159im, -0.0126871+0.0996538im, -0.269509+0.00302196im, -0.140838+0.0244051im, 0.000861337+0.00207717im, -6.42551e-5-0.00173344im, -0.00105528+0.00179777im, 0.000764672-0.000602481im, 0.00189714+0.00127471im, -0.30991+0.121537im, -0.0319424-0.144377im, -0.231995-0.0529414im, 0.241863-0.170549im, 0.0820029+0.0758833im, 0.168514+0.0935542im, 0.117687+0.0982952im, 0.0270232+0.081966im, -0.202623+0.122908im, -0.104325+0.0771085im, -0.156916-0.161661im, -0.101529-0.147743im, -0.000801824-0.117147im, 0.265459-0.0682365im, 0.135124-0.0486043im, 0.200206+0.268495im, 0.104333+0.207532im, -0.00332133+0.137411im, -0.382703+0.0442051im, -0.202173+0.0530631im, 0.00151036+0.00282271im, -0.000333829-0.00244654im, -0.00124394+0.002695im, 0.00099978-0.000960587im, 0.00286522+0.00154189im, 0.0811625-0.136874im, 0.286116+0.0340031im, 0.039079-0.115171im, -0.0396969+0.373233im, -0.160445+0.225487im, 0.12465+0.00265451im, 0.0853333+0.0120395im, 0.0583849+0.0470931im, -0.0738299+0.118299im, -0.0280568+0.0600721im, -0.386529+0.140124im, -0.270394+0.0470919im, -0.154767-0.0447004im, 0.0837792-0.458016im, 0.0072124-0.2535im, -0.133538+0.0105836im, -0.084408-0.00812938im, -0.0475752-0.0328252im, 0.0654166-0.130799im, 0.0203034-0.0724504im, -0.00122441-9.18638e-5im, 0.000860446+0.000395458im, -0.000616071-0.000956551im, 0.000105274+0.000520327im, -0.00108792+0.000611516im, 0.0298865-0.142912im, -0.217501+0.180601im, 0.175454-0.030206im, -0.13986-0.160109im, -0.0507848-0.328219im, 0.0525566-0.337304im, 0.0799062-0.22303im, 0.104006-0.0939662im, 0.304141+0.249739im, 0.180722+0.109822im, -0.0592788-0.0217604im, -0.0322979-0.0333287im, 0.013151-0.0260831im, 0.0586576-0.0412738im, 0.01969-0.0305178im, -0.134283+0.272541im, -0.121871+0.168365im, -0.113483+0.0536486im, -0.192178-0.28808im, -0.126192-0.138043im, -0.00164022+0.00237254im, 0.00182384-0.00127475im, -0.00267205+0.000112583im, 0.00118167+0.000405385im, -5.59629e-5+0.00293166im, -0.0899756+0.542342im, 0.284935+0.0481453im, -0.359211+0.11573im, -0.188399+0.0301327im, -0.201842+0.371847im, 0.061353-0.172194im, 0.0567285-0.113556im, 0.0525202-0.0403404im, 0.142324+0.156115im, 0.0951062+0.0784393im, 0.122461-0.0716501im, 0.0975112-0.0334494im, 0.0685959+0.0160209im, -0.00270279+0.165052im, 0.0093008+0.0828786im, -0.141691+0.0383132im, -0.101863+0.00900411im, -0.0553673-0.0233559im, 0.0437576-0.163564im, 0.0112594-0.0913605im, -0.00139728+0.000160781im, 0.00105287+0.000259349im, -0.0008995-0.000943779im, 0.000230458+0.000563954im, -0.00109281+0.000920657im], (ℂ^5 ⊗ ℂ^5) ← ℂ^5) +│ ⋮ + +```` + +Now that we have the state, we may compute observables using the `expectation_value` function. It typically expects a `Pair`, `(i1, i2, .., ik) => op` where `op` is a `TensorMap` or `InfiniteMPO` acting over `k` sites. In case of the Hamiltonian, it is not necessary to specify the indices as it spans the whole lattice. We can now plot the correlation function $\langle \hat{a}^{\dagger}_i \hat{a}_j\rangle$. + +````julia +plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:50), lw=2, xlabel="Site index", ylabel="Correlation function", yscale=:log10) +hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls=:dash, c=:black) +```` + +```@raw html + +``` + +We see that the correlations drop off exponentially, indicating the existence of a gapped Mott insulating phase. Let us now shift our parameters to probe other phases. + +````julia +ground_state = get_ground_state(0.5, 0.2, cutoff, D; tol=1e-6, verbosity=2, maxiter=500) + +plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:100), lw=2, xlabel="Site index", ylabel="Correlation function", yscale=:log10, xscale=:log10) +hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls=:dash, c=:black) +```` + +```@raw html + +``` + +In this case, the correlation function drops off algebraically and eventually saturates as $\lim_{i \to \infty}\langle\hat{a}_i^{\dagger} \hat{a}_j\rangle ≈ \langle \hat{a}_i^{\dagger}\rangle \langle \hat{a}_j \rangle = |\langle a_i \rangle|^2 \neq 0$. This is a signature of long-range order and suggests the existence of a Bose-Einstein condensate. However, this is a bit odd since at zero temperature, the Bose Hubbard model is not expected to break any continuous symmetries ($U(1)$ in this case, corresponding to particle number conservation) due to the Mermin-Wagner theorem. The source of this contradiction lies in the fact that the true 1D superfluid ground state is an extended critical phase exhibiting algebraic decay, however, a finite bond-dimension MPS can only capture exponentially decaying correlations. As a result, the finite bond dimension effectively introduces a length scale into the system in a similar manner as finite-size effects. We can see this clearly by increasing the bond dimension. We also see that the correlation length seems to depend algebraically on the bond dimension as expected from finite-entanglement scaling arguments. + +````julia +cutoff = 4 +Ds = 20:5:50 +mu, t = 0.5, 0.2 +states = Vector{InfiniteMPS}(undef, length(Ds)) + +Threads.@threads for idx in eachindex(Ds) + states[idx] = get_ground_state(mu, t, cutoff, Ds[idx], tol=1e-7, verbosity=1, maxiter=500) +end + +npoints = 400 +two_point_correlation = zeros(length(Ds), npoints) +a_op = a_min(cutoff=cutoff) + +Threads.@threads for idx in eachindex(Ds) + two_point_correlation[idx, :] .= real.(expectation_value(states[idx], (1, i) => a_op' ⊗ a_op) for i in 1:npoints) +end + +p = plot(framestyle=:box, ylabel="Correlation function " * L"\langle a_i^{\dagger}a_j \rangle", + xlabel="Distance " * L"|i-j|", xscale=:log10, yscale=:log10, + xticks=([10, 100], ["10", "100"]), + yticks=([0.25, 0.5, 1.0], ["0.25", "0.5", "1.0"])) + +plot!(p, 2:npoints, two_point_correlation[:, 2:end]', + lab="D = " .* string.(permutedims(Ds)), lw=2) + +scatter!(p, Ds, correlation_length.(states), + ylabel="Correlation length", xlabel="Bond dimension", + xscale=:log10, yscale=:log10, + inset=bbox(0.2, 0.51, 0.25, 0.25), + subplot=2, + xticks=(20:10:50, string.(20:10:50)), + yticks=([50, 100], string.([50, 100])), + xlabelfontsize=8, + ylabelfontsize=8, + ylims=[20, 130], + xlims=[15, 60]) +```` + +```@raw html + +``` + +This shows that any finite bond dimension MPS necessarily breaks the symmetry of the system, forming a Bose-Einstein condensate which introduces erraneous long-distance behaviour of correlation functions. In case of finite bond dimension, it is thus reasonable to associate the finite expectation value of the field operator to the 'quasicondensate' density of the system which vanishes as $D \to \infty$. + +````julia +quasicondensate_density = map(state -> abs2(expectation_value(state, (0,) => a_op)), states) +```` + +```` +7-element Vector{Float64}: + 0.3097427784562408 + 0.2881477510455394 + 0.2702002091214058 + 0.25712728499126414 + 0.24685384625099716 + 0.23539796049808664 + 0.22799690146561793 +```` + +We may now also visualize the momentum distribution function, which is obtained as the Fourier transform of the single-particle density matrix. Starting from the definition of the momentum occupation operator: + +$$\hat{a}_k = \frac{1}{\sqrt{L}} \sum_j e^{-ikj} \hat{a}_j, \qquad +\hat{a}_k^\dagger = \frac{1}{\sqrt{L}} \sum_{j'} e^{ikj'} \hat{a}_{j'}^\dagger +$$ + +the momentum distribution is + +$$\langle \hat{n}_k \rangle = \langle \hat{a}_k^\dagger \hat{a}_k \rangle += \frac{1}{L} \sum_{j',j} e^{ik(j'-j)} \langle \hat{a}_{j'}^\dagger \hat{a}_j \rangle.$$ + +For a translationally invariant system, the correlation depends only on the distance $r = j' - j$: + +$$\langle \hat{a}_{j'}^\dagger \hat{a}_j \rangle = C(r) = \langle \hat{a}_r^\dagger \hat{a}_0 \rangle.$$ + +Changing variables ($j' = j + r$) gives + +$$\langle \hat{n}_k \rangle = \frac{1}{L} \sum_j \sum_r e^{ikr} C(r).$$ + +The sum over $j$ yields a factor of $L$, which cancels the prefactor, leading to + +$$\langle \hat{n}_k \rangle = \sum_{r \in \mathbb{Z}} e^{ikr} \langle \hat{a}_r^\dagger \hat{a}_0 \rangle$$ + +However, we know that a finite bond dimension MPS introduces a non-zero quasi-condensate density which would give rise to an $\mathcal{O}(N)$ divergence in the momentum distribution that is not indicative of the true physics of the system. Since we know this contribution vanishes in the infinite bond dimension limit, we instead work with $\langle \hat{a}_r^{\dagger} \hat{a}_0 \rangle_c = \langle \hat{a}_r^{\dagger} \hat{a}_0 \rangle - |\langle \hat{a}\rangle|^2$. + +````julia +ks = range(-0.05, 0.15, 500) +momentum_distribution = map( + ((corr, qc),) -> sum( + 2 .* cos.(ks' .* (2:npoints)) .* (corr[2:end] .- qc), dims=1) .+ (corr[1] .- qc), + zip(eachrow(two_point_correlation), + quasicondensate_density)) +momentum_distribution = vcat(momentum_distribution...)' +plot(ks, momentum_distribution, lab="D = " .* string.(permutedims(Ds)), lw=1.5, xlabel="Momentum k", ylabel=L"\langle n_k \rangle", ylim=[0, 50]) +```` + +```@raw html + +``` + +We see that the density seems to peak around $k=0$, this time seemingly becoming more prominent as $D \to \infty$ which seems to suggest again that there is a condensate. However, going by the Penrose-Onsager criterion, the existence of a condensate can be quantified by requiring the leading eigenvalue of the single particle density matrix (i.e, $\langle \hat{n}_{k=0}\rangle = \sum_j \langle \hat{a}_j^{\dagger} \hat{a}_0\rangle$) to diverge as $O(N)$ in the thermodynamic limit. In this case, since the correlations decay as a power law, there is naturally a divergence at low momenta. But this does not imply the existence of a condensate since the order of divergence is much weaker. However, this does indicate the remnants of some kind of condensation in the 1D model despite the quantum fluctuations, leading to the practical utility of defining the concept of a quasicondensate where there is still a notion of phase coherence over short distances. + +What this means for us is that, as far as MPS simulations go, we may still utilize the quasicondensate density as an effective order parameter, although it will be less robust as the bond dimension is increased. Alternatively, we realize that the true phase is characterized as being a superfluid (a concept distinct from Bose-Einstein condensation) and can be identified by a non-zero value of the superfluid stiffness (also known as helicity modulus, $\Upsilon$) as defined by Leggett. Upon applying a phase twist $\Phi$ to the boundaries of the system, a superfluid phase would suffer an increase in energy whereas an insulating phase would not. In the thermodynamic limit, one could show that the boundary conditions may be considered as periodic and instead uniformly distribute the phase across the chain as $\hat{a}_i \to \hat{a}_i e^{i\Phi/L}$. Concretely, in the limit of $\Phi/L \to 0$, we have: + +$$\frac{E[\Phi] - E[0]}{L} \approx \frac{1}{2} \Upsilon(L) \bigg (\frac{\Phi}{L}\bigg)^2 + \cdots$$ + +In order to find the ground state under these twisted boundary conditions, we must construct our own variant of the Bose-Hubbard Hamiltonian. Typically you would want to take a peek at the [source code](https://github.com/QuantumKitHub/MPSKitModels.jl/blob/main/src/models/hamiltonians.jl#L379) of `MPSKitModels.jl` to see how these models are defined and tweak it as per your needs. Here we see that applying twisted boundary conditions is equivalent to adding a prefactor of $e^{\pm i\phi}$ in front of the hopping amplitudes. + +````julia +function bose_hubbard_model_twisted_bc(elt::Type{<:Number}=ComplexF64, symmetry::Type{<:Sector}=Trivial, + lattice::AbstractLattice=InfiniteChain(1); + cutoff::Integer=5, t=1.0, U=1.0, mu=0.0, phi=0) + + a_pm = a_plusmin(elt, symmetry; cutoff=cutoff) + a_mp = a_minplus(elt, symmetry; cutoff=cutoff) + N = a_number(elt, symmetry; cutoff=cutoff) + + interaction_term = N * (N - id(domain(N))) + + H = @mpoham begin + sum(nearest_neighbours(lattice)) do (i, j) + return -t * (exp(1im * phi) * a_pm{i,j} + exp(1im * -phi) * a_mp{i,j}) + end + + sum(vertices(lattice)) do i + return U / 2 * interaction_term{i} - mu * N{i} + end + end +end + +function superfluid_stiffness_profile(t, mu, D, cutoff, ϵ=1e-4, npoints=11) + phis = range(-ϵ, ϵ, npoints) + energies = zeros(length(phis)) + + Threads.@threads for idx in eachindex(phis) + hamiltonian_twisted = bose_hubbard_model_twisted_bc(cutoff=cutoff, t=t, mu=mu, U=1, phi=phis[idx]) + state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) + state_twisted, _, _ = find_groundstate(state, hamiltonian_twisted, VUMPS(; tol=1e-8, verbosity=0)) + energies[idx] = real(expectation_value(state_twisted, hamiltonian_twisted)) + end + + plot(phis, energies, lw=2, xlabel="Phase twist per site" * L"(\phi)", ylabel="Ground state energy", title="t = $t | μ = $mu | D = $D | cutoff = $cutoff") +end + +superfluid_stiffness_profile(0.2, 0.3, 5, 4) # superfluid + +superfluid_stiffness_profile(0.01, 0.3, 5, 4) # mott insulator +```` + +```@raw html + +``` + +Now that we know what phases to expect, we can plot the phase diagram by scanning over a range of parameters. In general, one could do better by performing a bisection algorithm for each chemical potential to determine the value of the hopping parameter at the transition point, however the 1D Bose-Hubbard model may have two transition points at the same chemical potential which makes this a bit cumbersome to implement robustly. Furthermore, we stick to using the quasi-condensate density as an order parameter since extracting the superfluid density accurately requires a more robust scheme to compute second derivatives which takes us away from the focus of this tutorial. + +````julia +cutoff, D = 4, 10 +mus = range(0, 0.75, 40) +ts = range(0, 0.3, 40) + +a_op = a_min(cutoff=cutoff) +order_parameters = zeros(length(ts), length(mus)) + +Threads.@threads for (i, j) in collect(Iterators.product(eachindex(mus), eachindex(ts))) + hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff=cutoff, U=1, mu=mus[i], t=ts[j]) + init_state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) + state, _, _ = find_groundstate(init_state, hamiltonian, VUMPS(; tol=1e-8, verbosity=0)) + order_parameters[i, j] = abs(expectation_value(state, 0 => a_op)) +end + +heatmap(ts, mus, order_parameters, xlabel=L"t/U", ylabel=L"\mu/U", title=L"\langle \hat{a}_i \rangle") +```` + +```@raw html + +``` + +Although the bond dimension here is quite low, we already see the deformation of the Mott insulator lobes to give way to the well known BKT transition that happens at commensurate density. One can go further and estimate the critical exponents using finite-entanglement scaling procedures on the correlation functions, but these may now be performed with ease using what we have learnt in this tutorial. + +--- + +*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* + diff --git a/docs/src/examples/quantum1d/8.bose-hubbard/main.ipynb b/docs/src/examples/quantum1d/8.bose-hubbard/main.ipynb new file mode 100644 index 000000000..b58889206 --- /dev/null +++ b/docs/src/examples/quantum1d/8.bose-hubbard/main.ipynb @@ -0,0 +1,424 @@ +{ + "cells": [ + { + "outputs": [], + "cell_type": "code", + "source": [ + "using Markdown\n", + "using MPSKit, MPSKitModels, TensorKit\n", + "using Plots, LaTeXStrings\n", + "\n", + "theme(:wong)\n", + "default(fontfamily=\"Computer Modern\", label=nothing, dpi=100, framestyle=:box)" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "# 1D Bose-Hubbard model\n", + "\n", + "In this tutorial, we will explore the physics of the one-dimensional Bose–Hubbard model using matrix product states. For the most part, we replicate the results presented in [**Phys. Rev. B 105, 134502**](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.105.134502), which can be consulted for any statements in this tutorial that are not otherwise cited. The Hamiltonian under study is now defined as follows:\n", + "\n", + "$$H = -t \\sum_{i} (\\hat{a}_i^{\\dagger} \\hat{a}_{i+1} + \\hat{a}_{i+1}^{\\dagger} \\hat{a}_i) + \\frac{U}{2} \\sum_i \\hat{n}_i(\\hat{n}_i - 1) - \\mu \\sum_i \\hat{n}_i$$\n", + "\n", + "where the bosonic creation and annihilation operators satisfy the canonical commutation relations (CCR):\n", + "\n", + "$$[\\hat{a}_i, \\hat{a}_j^{\\dagger}] = \\delta_{ij}.$$\n", + "\n", + "Each lattice site hosts a local Hilbert space corresponding to bosonic occupation states $|n\\rangle$, where $(n = 0, 1, 2, \\ldots)$. Since this space is formally infinite-dimensional, numerical simulations typically impose a truncation at some maximum occupation number $(n_{\\text{max}})$. Such a treatment is justified since it can be observed that the simulation results quickly converge with the cutoff if the filling fraction is kept sufficiently low.\n", + "\n", + "Within this truncated space, the local creation and annihilation operators are represented by finite-dimensional matrices. For example, with cutoff $n_{\\text{max}}$, the annihilation operator takes the form\n", + "\n", + "$$\n", + "\\hat{a} =\n", + "\\begin{bmatrix}\n", + "0 & \\sqrt{1} & 0 & 0 & \\cdots & 0 \\\\\n", + "0 & 0 & \\sqrt{2} & 0 & \\cdots & 0 \\\\\n", + "0 & 0 & 0 & \\sqrt{3} & \\cdots & 0 \\\\\n", + "\\vdots & & & \\ddots & \\ddots & \\vdots \\\\\n", + "0 & 0 & 0 & \\cdots & 0 & \\sqrt{n_{\\text{max}}} \\\\\n", + "0 & 0 & 0 & \\cdots & 0 & 0\n", + "\\end{bmatrix}\n", + "$$\n", + "\n", + "and the creation operator is simply its Hermitian conjugate,\n", + "\n", + "$$\n", + "\\hat{a}^\\dagger =\n", + "\\begin{bmatrix}\n", + "0 & 0 & 0 & \\cdots & 0 & 0 \\\\\n", + "\\sqrt{1} & 0 & 0 & \\cdots & 0 & 0 \\\\\n", + "0 & \\sqrt{2} & 0 & \\cdots & 0 & 0 \\\\\n", + "\\vdots & & \\ddots & \\ddots & & \\vdots \\\\\n", + "0 & 0 & 0 & \\cdots & 0 & 0 \\\\\n", + "0 & 0 & 0 & \\cdots & \\sqrt{n_{\\text{max}}} & 0\n", + "\\end{bmatrix}\n", + "$$\n", + "\n", + "The number operator is then given by\n", + "\n", + "$$\\hat{n} = \\hat{a}^\\dagger \\hat{a} = \\mathrm{diag}(0, 1, 2, \\ldots, n_{\\text{max}}).$$\n", + "\n", + "Before moving on, notice that the Hamiltonian is uniform and translationally invariant. Typically, such models are studied on a finite chain of $N$ sites with periodic boundary conditions, but this introduces finite-size effects that are rather annoying to deal with. In contrast, the MPS framework allows us to work directly in the thermodynamic limit, avoiding such artifacts. We will follow this line of exploration in this tutorial and leave finite systems for another time.\n", + "\n", + "In order to work in the thermodynamic limit, we will have to create an `InfiniteMPS`. A complete specification of the MPS requires us to define the physical space and the virtual space of the constituent tensors. At this point is it useful to note that `MPSKit.jl` is powered by `TensorKit.jl` under the hood and has some very generic interfaces in order to allow imposing symmetries of all kinds. As a result, it is sometimes necessary to be a bit more explicit about what we want to do in terms of the vector spaces involved. In this case, we will not consider any symmetries and simply take the most naive approach of working within the `Trivial` sector. The physical space is then `ℂ^(nmax+1)`, and the virtual space is `ℂ^D` where $D$ is some integer chosen to be the bond dimension of the MPS, and `ℂ` is an alias for `ComplexSpace` (typeset as `\\bbC`). As $D$ is increased, one increases the amount of entanglement, i.e, quantum correlations that can be captured by the state." + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "cutoff, D = 4, 5\n", + "initial_state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D)" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "This simply initializes a tensor filled with random entries (check out the documentation for other useful constructors). Next, we need the creation and annihilation operators. While we could construct them from scratch, here we will use `MPSKitModels.jl` instead that has predefined operators and models for most well-known lattice models." + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "a_op = a_min(cutoff=cutoff) # creates a bosonic annihilation operator without any symmetries\n", + "display(a_op[])\n", + "display((a_op'*a_op)[])" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "The [] accessor lets us see the underlying array and indeed the operators are exactly what we require. Similarly, the Bose Hubbard model is also predefined (although we will construct our own variant later on)." + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "hamiltonian = bose_hubbard_model(InfiniteChain(1), cutoff=cutoff, U=1, mu=0.5, t=0.2) # It is not strictly required to pass InfiniteChain() and is only included for clarity; one may instead pass FiniteChain(N) as well" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "This has created the Hamiltonian operator as a matrix product operator (MPO) which is a convenient form to use in conjunction with MPS. Finally, the ground state optimization may be performed with either `iDMRG` or `VUMPS`. Both should take similar arguments but it is known that VUMPS is typically more efficient for these systems so we proceed with that." + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "ground_state, _, _ = find_groundstate(initial_state, hamiltonian, VUMPS(tol=1e-6, verbosity=2, maxiter=200))\n", + "println(\"Energy: \", expectation_value(ground_state, hamiltonian))" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "This automatically runs the algorithm until a certain error [measure](https://github.com/QuantumKitHub/MPSKit.jl/blob/main/src/algorithms/toolbox.jl#L42) falls below the specified tolerance or the maximum iterations is reached. Let us wrap all this into a convenient function." + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "function get_ground_state(mu, t, cutoff, D; kwargs...)\n", + " hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff=cutoff, U=1, mu=mu, t=t)\n", + " state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D)\n", + " state, _, _ = find_groundstate(state, hamiltonian, VUMPS(; kwargs...))\n", + "\n", + " return state\n", + "end\n", + "\n", + "ground_state = get_ground_state(0.5, 0.01, cutoff, D; tol=1e-6, verbosity=2, maxiter=500)" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "Now that we have the state, we may compute observables using the `expectation_value` function. It typically expects a `Pair`, `(i1, i2, .., ik) => op` where `op` is a `TensorMap` or `InfiniteMPO` acting over `k` sites. In case of the Hamiltonian, it is not necessary to specify the indices as it spans the whole lattice. We can now plot the correlation function $\\langle \\hat{a}^{\\dagger}_i \\hat{a}_j\\rangle$." + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:50), lw=2, xlabel=\"Site index\", ylabel=\"Correlation function\", yscale=:log10)\n", + "hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls=:dash, c=:black)" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "We see that the correlations drop off exponentially, indicating the existence of a gapped Mott insulating phase. Let us now shift our parameters to probe other phases." + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "ground_state = get_ground_state(0.5, 0.2, cutoff, D; tol=1e-6, verbosity=2, maxiter=500)\n", + "\n", + "plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:100), lw=2, xlabel=\"Site index\", ylabel=\"Correlation function\", yscale=:log10, xscale=:log10)\n", + "hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls=:dash, c=:black)" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "In this case, the correlation function drops off algebraically and eventually saturates as $\\lim_{i \\to \\infty}\\langle\\hat{a}_i^{\\dagger} \\hat{a}_j\\rangle ≈ \\langle \\hat{a}_i^{\\dagger}\\rangle \\langle \\hat{a}_j \\rangle = |\\langle a_i \\rangle|^2 \\neq 0$. This is a signature of long-range order and suggests the existence of a Bose-Einstein condensate. However, this is a bit odd since at zero temperature, the Bose Hubbard model is not expected to break any continuous symmetries ($U(1)$ in this case, corresponding to particle number conservation) due to the Mermin-Wagner theorem. The source of this contradiction lies in the fact that the true 1D superfluid ground state is an extended critical phase exhibiting algebraic decay, however, a finite bond-dimension MPS can only capture exponentially decaying correlations. As a result, the finite bond dimension effectively introduces a length scale into the system in a similar manner as finite-size effects. We can see this clearly by increasing the bond dimension. We also see that the correlation length seems to depend algebraically on the bond dimension as expected from finite-entanglement scaling arguments." + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "cutoff = 4\n", + "Ds = 20:5:50\n", + "mu, t = 0.5, 0.2\n", + "states = Vector{InfiniteMPS}(undef, length(Ds))\n", + "\n", + "Threads.@threads for idx in eachindex(Ds)\n", + " states[idx] = get_ground_state(mu, t, cutoff, Ds[idx], tol=1e-7, verbosity=1, maxiter=500)\n", + "end\n", + "\n", + "npoints = 400\n", + "two_point_correlation = zeros(length(Ds), npoints)\n", + "a_op = a_min(cutoff=cutoff)\n", + "\n", + "Threads.@threads for idx in eachindex(Ds)\n", + " two_point_correlation[idx, :] .= real.(expectation_value(states[idx], (1, i) => a_op' ⊗ a_op) for i in 1:npoints)\n", + "end\n", + "\n", + "p = plot(framestyle=:box, ylabel=\"Correlation function \" * L\"\\langle a_i^{\\dagger}a_j \\rangle\",\n", + " xlabel=\"Distance \" * L\"|i-j|\", xscale=:log10, yscale=:log10,\n", + " xticks=([10, 100], [\"10\", \"100\"]),\n", + " yticks=([0.25, 0.5, 1.0], [\"0.25\", \"0.5\", \"1.0\"]))\n", + "\n", + "plot!(p, 2:npoints, two_point_correlation[:, 2:end]',\n", + " lab=\"D = \" .* string.(permutedims(Ds)), lw=2)\n", + "\n", + "scatter!(p, Ds, correlation_length.(states),\n", + " ylabel=\"Correlation length\", xlabel=\"Bond dimension\",\n", + " xscale=:log10, yscale=:log10,\n", + " inset=bbox(0.2, 0.51, 0.25, 0.25),\n", + " subplot=2,\n", + " xticks=(20:10:50, string.(20:10:50)),\n", + " yticks=([50, 100], string.([50, 100])),\n", + " xlabelfontsize=8,\n", + " ylabelfontsize=8,\n", + " ylims=[20, 130],\n", + " xlims=[15, 60])" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "This shows that any finite bond dimension MPS necessarily breaks the symmetry of the system, forming a Bose-Einstein condensate which introduces erraneous long-distance behaviour of correlation functions. In case of finite bond dimension, it is thus reasonable to associate the finite expectation value of the field operator to the 'quasicondensate' density of the system which vanishes as $D \\to \\infty$." + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "quasicondensate_density = map(state -> abs2(expectation_value(state, (0,) => a_op)), states)" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "We may now also visualize the momentum distribution function, which is obtained as the Fourier transform of the single-particle density matrix. Starting from the definition of the momentum occupation operator:\n", + "\n", + "$$\\hat{a}_k = \\frac{1}{\\sqrt{L}} \\sum_j e^{-ikj} \\hat{a}_j, \\qquad\n", + "\\hat{a}_k^\\dagger = \\frac{1}{\\sqrt{L}} \\sum_{j'} e^{ikj'} \\hat{a}_{j'}^\\dagger\n", + "$$\n", + "\n", + "the momentum distribution is\n", + "\n", + "$$\\langle \\hat{n}_k \\rangle = \\langle \\hat{a}_k^\\dagger \\hat{a}_k \\rangle\n", + "= \\frac{1}{L} \\sum_{j',j} e^{ik(j'-j)} \\langle \\hat{a}_{j'}^\\dagger \\hat{a}_j \\rangle.$$\n", + "\n", + "For a translationally invariant system, the correlation depends only on the distance $r = j' - j$:\n", + "\n", + "$$\\langle \\hat{a}_{j'}^\\dagger \\hat{a}_j \\rangle = C(r) = \\langle \\hat{a}_r^\\dagger \\hat{a}_0 \\rangle.$$\n", + "\n", + "Changing variables ($j' = j + r$) gives\n", + "\n", + "$$\\langle \\hat{n}_k \\rangle = \\frac{1}{L} \\sum_j \\sum_r e^{ikr} C(r).$$\n", + "\n", + "The sum over $j$ yields a factor of $L$, which cancels the prefactor, leading to\n", + "\n", + "$$\\langle \\hat{n}_k \\rangle = \\sum_{r \\in \\mathbb{Z}} e^{ikr} \\langle \\hat{a}_r^\\dagger \\hat{a}_0 \\rangle$$\n", + "\n", + "However, we know that a finite bond dimension MPS introduces a non-zero quasi-condensate density which would give rise to an $\\mathcal{O}(N)$ divergence in the momentum distribution that is not indicative of the true physics of the system. Since we know this contribution vanishes in the infinite bond dimension limit, we instead work with $\\langle \\hat{a}_r^{\\dagger} \\hat{a}_0 \\rangle_c = \\langle \\hat{a}_r^{\\dagger} \\hat{a}_0 \\rangle - |\\langle \\hat{a}\\rangle|^2$." + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "ks = range(-0.05, 0.15, 500)\n", + "momentum_distribution = map(\n", + " ((corr, qc),) -> sum(\n", + " 2 .* cos.(ks' .* (2:npoints)) .* (corr[2:end] .- qc), dims=1) .+ (corr[1] .- qc),\n", + " zip(eachrow(two_point_correlation),\n", + " quasicondensate_density))\n", + "momentum_distribution = vcat(momentum_distribution...)'\n", + "plot(ks, momentum_distribution, lab=\"D = \" .* string.(permutedims(Ds)), lw=1.5, xlabel=\"Momentum k\", ylabel=L\"\\langle n_k \\rangle\", ylim=[0, 50])" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "We see that the density seems to peak around $k=0$, this time seemingly becoming more prominent as $D \\to \\infty$ which seems to suggest again that there is a condensate. However, going by the Penrose-Onsager criterion, the existence of a condensate can be quantified by requiring the leading eigenvalue of the single particle density matrix (i.e, $\\langle \\hat{n}_{k=0}\\rangle = \\sum_j \\langle \\hat{a}_j^{\\dagger} \\hat{a}_0\\rangle$) to diverge as $O(N)$ in the thermodynamic limit. In this case, since the correlations decay as a power law, there is naturally a divergence at low momenta. But this does not imply the existence of a condensate since the order of divergence is much weaker. However, this does indicate the remnants of some kind of condensation in the 1D model despite the quantum fluctuations, leading to the practical utility of defining the concept of a quasicondensate where there is still a notion of phase coherence over short distances.\n", + "\n", + "What this means for us is that, as far as MPS simulations go, we may still utilize the quasicondensate density as an effective order parameter, although it will be less robust as the bond dimension is increased. Alternatively, we realize that the true phase is characterized as being a superfluid (a concept distinct from Bose-Einstein condensation) and can be identified by a non-zero value of the superfluid stiffness (also known as helicity modulus, $\\Upsilon$) as defined by Leggett. Upon applying a phase twist $\\Phi$ to the boundaries of the system, a superfluid phase would suffer an increase in energy whereas an insulating phase would not. In the thermodynamic limit, one could show that the boundary conditions may be considered as periodic and instead uniformly distribute the phase across the chain as $\\hat{a}_i \\to \\hat{a}_i e^{i\\Phi/L}$. Concretely, in the limit of $\\Phi/L \\to 0$, we have:\n", + "\n", + "$$\\frac{E[\\Phi] - E[0]}{L} \\approx \\frac{1}{2} \\Upsilon(L) \\bigg (\\frac{\\Phi}{L}\\bigg)^2 + \\cdots$$\n", + "\n", + "In order to find the ground state under these twisted boundary conditions, we must construct our own variant of the Bose-Hubbard Hamiltonian. Typically you would want to take a peek at the [source code](https://github.com/QuantumKitHub/MPSKitModels.jl/blob/main/src/models/hamiltonians.jl#L379) of `MPSKitModels.jl` to see how these models are defined and tweak it as per your needs. Here we see that applying twisted boundary conditions is equivalent to adding a prefactor of $e^{\\pm i\\phi}$ in front of the hopping amplitudes." + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "function bose_hubbard_model_twisted_bc(elt::Type{<:Number}=ComplexF64, symmetry::Type{<:Sector}=Trivial,\n", + " lattice::AbstractLattice=InfiniteChain(1);\n", + " cutoff::Integer=5, t=1.0, U=1.0, mu=0.0, phi=0)\n", + "\n", + " a_pm = a_plusmin(elt, symmetry; cutoff=cutoff)\n", + " a_mp = a_minplus(elt, symmetry; cutoff=cutoff)\n", + " N = a_number(elt, symmetry; cutoff=cutoff)\n", + "\n", + " interaction_term = N * (N - id(domain(N)))\n", + "\n", + " H = @mpoham begin\n", + " sum(nearest_neighbours(lattice)) do (i, j)\n", + " return -t * (exp(1im * phi) * a_pm{i,j} + exp(1im * -phi) * a_mp{i,j})\n", + " end +\n", + " sum(vertices(lattice)) do i\n", + " return U / 2 * interaction_term{i} - mu * N{i}\n", + " end\n", + " end\n", + "end\n", + "\n", + "function superfluid_stiffness_profile(t, mu, D, cutoff, ϵ=1e-4, npoints=11)\n", + " phis = range(-ϵ, ϵ, npoints)\n", + " energies = zeros(length(phis))\n", + "\n", + " Threads.@threads for idx in eachindex(phis)\n", + " hamiltonian_twisted = bose_hubbard_model_twisted_bc(cutoff=cutoff, t=t, mu=mu, U=1, phi=phis[idx])\n", + " state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D)\n", + " state_twisted, _, _ = find_groundstate(state, hamiltonian_twisted, VUMPS(; tol=1e-8, verbosity=0))\n", + " energies[idx] = real(expectation_value(state_twisted, hamiltonian_twisted))\n", + " end\n", + "\n", + " plot(phis, energies, lw=2, xlabel=\"Phase twist per site\" * L\"(\\phi)\", ylabel=\"Ground state energy\", title=\"t = $t | μ = $mu | D = $D | cutoff = $cutoff\")\n", + "end\n", + "\n", + "superfluid_stiffness_profile(0.2, 0.3, 5, 4) # superfluid\n", + "\n", + "superfluid_stiffness_profile(0.01, 0.3, 5, 4) # mott insulator" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "Now that we know what phases to expect, we can plot the phase diagram by scanning over a range of parameters. In general, one could do better by performing a bisection algorithm for each chemical potential to determine the value of the hopping parameter at the transition point, however the 1D Bose-Hubbard model may have two transition points at the same chemical potential which makes this a bit cumbersome to implement robustly. Furthermore, we stick to using the quasi-condensate density as an order parameter since extracting the superfluid density accurately requires a more robust scheme to compute second derivatives which takes us away from the focus of this tutorial." + ], + "metadata": {} + }, + { + "outputs": [], + "cell_type": "code", + "source": [ + "cutoff, D = 4, 10\n", + "mus = range(0, 0.75, 40)\n", + "ts = range(0, 0.3, 40)\n", + "\n", + "a_op = a_min(cutoff=cutoff)\n", + "order_parameters = zeros(length(ts), length(mus))\n", + "\n", + "Threads.@threads for (i, j) in collect(Iterators.product(eachindex(mus), eachindex(ts)))\n", + " hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff=cutoff, U=1, mu=mus[i], t=ts[j])\n", + " init_state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D)\n", + " state, _, _ = find_groundstate(init_state, hamiltonian, VUMPS(; tol=1e-8, verbosity=0))\n", + " order_parameters[i, j] = abs(expectation_value(state, 0 => a_op))\n", + "end\n", + "\n", + "heatmap(ts, mus, order_parameters, xlabel=L\"t/U\", ylabel=L\"\\mu/U\", title=L\"\\langle \\hat{a}_i \\rangle\")" + ], + "metadata": {}, + "execution_count": null + }, + { + "cell_type": "markdown", + "source": [ + "Although the bond dimension here is quite low, we already see the deformation of the Mott insulator lobes to give way to the well known BKT transition that happens at commensurate density. One can go further and estimate the critical exponents using finite-entanglement scaling procedures on the correlation functions, but these may now be performed with ease using what we have learnt in this tutorial." + ], + "metadata": {} + }, + { + "cell_type": "markdown", + "source": [ + "---\n", + "\n", + "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*" + ], + "metadata": {} + } + ], + "nbformat_minor": 3, + "metadata": { + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.11.7" + }, + "kernelspec": { + "name": "julia-1.11", + "display_name": "Julia 1.11.7", + "language": "julia" + } + }, + "nbformat": 4 +} \ No newline at end of file diff --git a/examples/Cache.toml b/examples/Cache.toml index 8820b4277..a4a34c75a 100644 --- a/examples/Cache.toml +++ b/examples/Cache.toml @@ -5,6 +5,7 @@ "2.haldane" = "804433b690faa1ce268a430edb4185af77a38de7a9f53e097795688a53c30c5e" "6.hubbard" = "af14e1df11b2392a69260ce061a716370a2ce50b5c907f0a7d2548e796d543fe" "7.xy-finiteT" = "7afb26fb9ff8ef84722fee3442eaed2ddffda4a5f70a7844b61ce5178093706c" +"8.bose-hubbard" = "63c8df67cede6db14296a87aa86b2e62fde50d0ab5cde6401c07ca692d18ea99" "3.ising-dqpt" = "2d314fd05a75c5c91ff81391c7bf166171b35a7b32933e9490756dc584fe8437" "5.haldane-spt" = "3de1b3baa1e4c5dc2252bbe8688d0c6ac8e5d5265deeb6ab66818bbfb02dd4aa" "4.xxz-heisenberg" = "1a2093a182f3ce070b53488484d1024e45496b291c1b74e3f370201d4c056be2" diff --git a/examples/Project.toml b/examples/Project.toml index 5c2ff210a..3a4751664 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -2,6 +2,7 @@ BenchmarkFreeFermions = "5b68a00d-ca4d-4b58-ab0c-dc082ade624e" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" KrylovKit = "0b1a1467-8014-51b9-945f-bf0ae24f4b77" +LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" MPSKit = "bb1c41ca-d63c-52ed-829e-0820dda26502" diff --git a/examples/quantum1d/8.bose-hubbard/main.jl b/examples/quantum1d/8.bose-hubbard/main.jl new file mode 100644 index 000000000..dd7e2a3c8 --- /dev/null +++ b/examples/quantum1d/8.bose-hubbard/main.jl @@ -0,0 +1,265 @@ +using Markdown +using MPSKit, MPSKitModels, TensorKit +using Plots, LaTeXStrings + +theme(:wong) +default(fontfamily="Computer Modern", label=nothing, dpi=100, framestyle=:box) + +md""" +# 1D Bose-Hubbard model + +In this tutorial, we will explore the physics of the one-dimensional Bose–Hubbard model using matrix product states. For the most part, we replicate the results presented in [**Phys. Rev. B 105, 134502**](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.105.134502), which can be consulted for any statements in this tutorial that are not otherwise cited. The Hamiltonian under study is now defined as follows: + +$$H = -t \sum_{i} (\hat{a}_i^{\dagger} \hat{a}_{i+1} + \hat{a}_{i+1}^{\dagger} \hat{a}_i) + \frac{U}{2} \sum_i \hat{n}_i(\hat{n}_i - 1) - \mu \sum_i \hat{n}_i$$ + +where the bosonic creation and annihilation operators satisfy the canonical commutation relations (CCR): + +$$[\hat{a}_i, \hat{a}_j^{\dagger}] = \delta_{ij}.$$ + +Each lattice site hosts a local Hilbert space corresponding to bosonic occupation states $|n\rangle$, where $(n = 0, 1, 2, \ldots)$. Since this space is formally infinite-dimensional, numerical simulations typically impose a truncation at some maximum occupation number $(n_{\text{max}})$. Such a treatment is justified since it can be observed that the simulation results quickly converge with the cutoff if the filling fraction is kept sufficiently low. + +Within this truncated space, the local creation and annihilation operators are represented by finite-dimensional matrices. For example, with cutoff $n_{\text{max}}$, the annihilation operator takes the form + +$$ +\hat{a} = +\begin{bmatrix} +0 & \sqrt{1} & 0 & 0 & \cdots & 0 \\ +0 & 0 & \sqrt{2} & 0 & \cdots & 0 \\ +0 & 0 & 0 & \sqrt{3} & \cdots & 0 \\ +\vdots & & & \ddots & \ddots & \vdots \\ +0 & 0 & 0 & \cdots & 0 & \sqrt{n_{\text{max}}} \\ +0 & 0 & 0 & \cdots & 0 & 0 +\end{bmatrix} +$$ + +and the creation operator is simply its Hermitian conjugate, + +$$ +\hat{a}^\dagger = +\begin{bmatrix} +0 & 0 & 0 & \cdots & 0 & 0 \\ +\sqrt{1} & 0 & 0 & \cdots & 0 & 0 \\ +0 & \sqrt{2} & 0 & \cdots & 0 & 0 \\ +\vdots & & \ddots & \ddots & & \vdots \\ +0 & 0 & 0 & \cdots & 0 & 0 \\ +0 & 0 & 0 & \cdots & \sqrt{n_{\text{max}}} & 0 +\end{bmatrix} +$$ + +The number operator is then given by + +$$\hat{n} = \hat{a}^\dagger \hat{a} = \mathrm{diag}(0, 1, 2, \ldots, n_{\text{max}}).$$ + +Before moving on, notice that the Hamiltonian is uniform and translationally invariant. Typically, such models are studied on a finite chain of $N$ sites with periodic boundary conditions, but this introduces finite-size effects that are rather annoying to deal with. In contrast, the MPS framework allows us to work directly in the thermodynamic limit, avoiding such artifacts. We will follow this line of exploration in this tutorial and leave finite systems for another time. + +In order to work in the thermodynamic limit, we will have to create an `InfiniteMPS`. A complete specification of the MPS requires us to define the physical space and the virtual space of the constituent tensors. At this point is it useful to note that `MPSKit.jl` is powered by `TensorKit.jl` under the hood and has some very generic interfaces in order to allow imposing symmetries of all kinds. As a result, it is sometimes necessary to be a bit more explicit about what we want to do in terms of the vector spaces involved. In this case, we will not consider any symmetries and simply take the most naive approach of working within the `Trivial` sector. The physical space is then `ℂ^(nmax+1)`, and the virtual space is `ℂ^D` where $D$ is some integer chosen to be the bond dimension of the MPS, and `ℂ` is an alias for `ComplexSpace` (typeset as `\bbC`). As $D$ is increased, one increases the amount of entanglement, i.e, quantum correlations that can be captured by the state. +""" + +cutoff, D = 4, 5 +initial_state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) + +md""" +This simply initializes a tensor filled with random entries (check out the documentation for other useful constructors). Next, we need the creation and annihilation operators. While we could construct them from scratch, here we will use `MPSKitModels.jl` instead that has predefined operators and models for most well-known lattice models. +""" + +a_op = a_min(cutoff=cutoff) # creates a bosonic annihilation operator without any symmetries +display(a_op[]) +display((a_op'*a_op)[]) + +md""" +The [] accessor lets us see the underlying array and indeed the operators are exactly what we require. Similarly, the Bose Hubbard model is also predefined (although we will construct our own variant later on). +""" + +hamiltonian = bose_hubbard_model(InfiniteChain(1), cutoff=cutoff, U=1, mu=0.5, t=0.2) # It is not strictly required to pass InfiniteChain() and is only included for clarity; one may instead pass FiniteChain(N) as well + +md""" +This has created the Hamiltonian operator as a matrix product operator (MPO) which is a convenient form to use in conjunction with MPS. Finally, the ground state optimization may be performed with either `iDMRG` or `VUMPS`. Both should take similar arguments but it is known that VUMPS is typically more efficient for these systems so we proceed with that. +""" + +ground_state, _, _ = find_groundstate(initial_state, hamiltonian, VUMPS(tol=1e-6, verbosity=2, maxiter=200)) +println("Energy: ", expectation_value(ground_state, hamiltonian)) + +md""" +This automatically runs the algorithm until a certain error [measure](https://github.com/QuantumKitHub/MPSKit.jl/blob/main/src/algorithms/toolbox.jl#L42) falls below the specified tolerance or the maximum iterations is reached. Let us wrap all this into a convenient function. +""" + +function get_ground_state(mu, t, cutoff, D; kwargs...) + hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff=cutoff, U=1, mu=mu, t=t) + state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) + state, _, _ = find_groundstate(state, hamiltonian, VUMPS(; kwargs...)) + + return state +end + +ground_state = get_ground_state(0.5, 0.01, cutoff, D; tol=1e-6, verbosity=2, maxiter=500) + +md""" +Now that we have the state, we may compute observables using the `expectation_value` function. It typically expects a `Pair`, `(i1, i2, .., ik) => op` where `op` is a `TensorMap` or `InfiniteMPO` acting over `k` sites. In case of the Hamiltonian, it is not necessary to specify the indices as it spans the whole lattice. We can now plot the correlation function $\langle \hat{a}^{\dagger}_i \hat{a}_j\rangle$. +""" + +plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:50), lw=2, xlabel="Site index", ylabel="Correlation function", yscale=:log10) +hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls=:dash, c=:black) + +md""" +We see that the correlations drop off exponentially, indicating the existence of a gapped Mott insulating phase. Let us now shift our parameters to probe other phases. +""" + +ground_state = get_ground_state(0.5, 0.2, cutoff, D; tol=1e-6, verbosity=2, maxiter=500) + +plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:100), lw=2, xlabel="Site index", ylabel="Correlation function", yscale=:log10, xscale=:log10) +hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls=:dash, c=:black) + +md""" +In this case, the correlation function drops off algebraically and eventually saturates as $\lim_{i \to \infty}\langle\hat{a}_i^{\dagger} \hat{a}_j\rangle ≈ \langle \hat{a}_i^{\dagger}\rangle \langle \hat{a}_j \rangle = |\langle a_i \rangle|^2 \neq 0$. This is a signature of long-range order and suggests the existence of a Bose-Einstein condensate. However, this is a bit odd since at zero temperature, the Bose Hubbard model is not expected to break any continuous symmetries ($U(1)$ in this case, corresponding to particle number conservation) due to the Mermin-Wagner theorem. The source of this contradiction lies in the fact that the true 1D superfluid ground state is an extended critical phase exhibiting algebraic decay, however, a finite bond-dimension MPS can only capture exponentially decaying correlations. As a result, the finite bond dimension effectively introduces a length scale into the system in a similar manner as finite-size effects. We can see this clearly by increasing the bond dimension. We also see that the correlation length seems to depend algebraically on the bond dimension as expected from finite-entanglement scaling arguments. +""" + +cutoff = 4 +Ds = 20:5:50 +mu, t = 0.5, 0.2 +states = Vector{InfiniteMPS}(undef, length(Ds)) + +Threads.@threads for idx in eachindex(Ds) + states[idx] = get_ground_state(mu, t, cutoff, Ds[idx], tol=1e-7, verbosity=1, maxiter=500) +end + +npoints = 400 +two_point_correlation = zeros(length(Ds), npoints) +a_op = a_min(cutoff=cutoff) + +Threads.@threads for idx in eachindex(Ds) + two_point_correlation[idx, :] .= real.(expectation_value(states[idx], (1, i) => a_op' ⊗ a_op) for i in 1:npoints) +end + +p = plot(framestyle=:box, ylabel="Correlation function " * L"\langle a_i^{\dagger}a_j \rangle", + xlabel="Distance " * L"|i-j|", xscale=:log10, yscale=:log10, + xticks=([10, 100], ["10", "100"]), + yticks=([0.25, 0.5, 1.0], ["0.25", "0.5", "1.0"])) + +plot!(p, 2:npoints, two_point_correlation[:, 2:end]', + lab="D = " .* string.(permutedims(Ds)), lw=2) + +scatter!(p, Ds, correlation_length.(states), + ylabel="Correlation length", xlabel="Bond dimension", + xscale=:log10, yscale=:log10, + inset=bbox(0.2, 0.51, 0.25, 0.25), + subplot=2, + xticks=(20:10:50, string.(20:10:50)), + yticks=([50, 100], string.([50, 100])), + xlabelfontsize=8, + ylabelfontsize=8, + ylims=[20, 130], + xlims=[15, 60]) + +md""" +This shows that any finite bond dimension MPS necessarily breaks the symmetry of the system, forming a Bose-Einstein condensate which introduces erraneous long-distance behaviour of correlation functions. In case of finite bond dimension, it is thus reasonable to associate the finite expectation value of the field operator to the 'quasicondensate' density of the system which vanishes as $D \to \infty$. +""" + +quasicondensate_density = map(state -> abs2(expectation_value(state, (0,) => a_op)), states) + +md""" +We may now also visualize the momentum distribution function, which is obtained as the Fourier transform of the single-particle density matrix. Starting from the definition of the momentum occupation operator: + +$$\hat{a}_k = \frac{1}{\sqrt{L}} \sum_j e^{-ikj} \hat{a}_j, \qquad +\hat{a}_k^\dagger = \frac{1}{\sqrt{L}} \sum_{j'} e^{ikj'} \hat{a}_{j'}^\dagger +$$ + +the momentum distribution is + +$$\langle \hat{n}_k \rangle = \langle \hat{a}_k^\dagger \hat{a}_k \rangle += \frac{1}{L} \sum_{j',j} e^{ik(j'-j)} \langle \hat{a}_{j'}^\dagger \hat{a}_j \rangle.$$ + +For a translationally invariant system, the correlation depends only on the distance $r = j' - j$: + +$$\langle \hat{a}_{j'}^\dagger \hat{a}_j \rangle = C(r) = \langle \hat{a}_r^\dagger \hat{a}_0 \rangle.$$ + +Changing variables ($j' = j + r$) gives + +$$\langle \hat{n}_k \rangle = \frac{1}{L} \sum_j \sum_r e^{ikr} C(r).$$ + +The sum over $j$ yields a factor of $L$, which cancels the prefactor, leading to + +$$\langle \hat{n}_k \rangle = \sum_{r \in \mathbb{Z}} e^{ikr} \langle \hat{a}_r^\dagger \hat{a}_0 \rangle$$ + +However, we know that a finite bond dimension MPS introduces a non-zero quasi-condensate density which would give rise to an $\mathcal{O}(N)$ divergence in the momentum distribution that is not indicative of the true physics of the system. Since we know this contribution vanishes in the infinite bond dimension limit, we instead work with $\langle \hat{a}_r^{\dagger} \hat{a}_0 \rangle_c = \langle \hat{a}_r^{\dagger} \hat{a}_0 \rangle - |\langle \hat{a}\rangle|^2$. +""" + +ks = range(-0.05, 0.15, 500) +momentum_distribution = map( + ((corr, qc),) -> sum( + 2 .* cos.(ks' .* (2:npoints)) .* (corr[2:end] .- qc), dims=1) .+ (corr[1] .- qc), + zip(eachrow(two_point_correlation), + quasicondensate_density)) +momentum_distribution = vcat(momentum_distribution...)' +plot(ks, momentum_distribution, lab="D = " .* string.(permutedims(Ds)), lw=1.5, xlabel="Momentum k", ylabel=L"\langle n_k \rangle", ylim=[0, 50]) + +md""" +We see that the density seems to peak around $k=0$, this time seemingly becoming more prominent as $D \to \infty$ which seems to suggest again that there is a condensate. However, going by the Penrose-Onsager criterion, the existence of a condensate can be quantified by requiring the leading eigenvalue of the single particle density matrix (i.e, $\langle \hat{n}_{k=0}\rangle = \sum_j \langle \hat{a}_j^{\dagger} \hat{a}_0\rangle$) to diverge as $O(N)$ in the thermodynamic limit. In this case, since the correlations decay as a power law, there is naturally a divergence at low momenta. But this does not imply the existence of a condensate since the order of divergence is much weaker. However, this does indicate the remnants of some kind of condensation in the 1D model despite the quantum fluctuations, leading to the practical utility of defining the concept of a quasicondensate where there is still a notion of phase coherence over short distances. + +What this means for us is that, as far as MPS simulations go, we may still utilize the quasicondensate density as an effective order parameter, although it will be less robust as the bond dimension is increased. Alternatively, we realize that the true phase is characterized as being a superfluid (a concept distinct from Bose-Einstein condensation) and can be identified by a non-zero value of the superfluid stiffness (also known as helicity modulus, $\Upsilon$) as defined by Leggett. Upon applying a phase twist $\Phi$ to the boundaries of the system, a superfluid phase would suffer an increase in energy whereas an insulating phase would not. In the thermodynamic limit, one could show that the boundary conditions may be considered as periodic and instead uniformly distribute the phase across the chain as $\hat{a}_i \to \hat{a}_i e^{i\Phi/L}$. Concretely, in the limit of $\Phi/L \to 0$, we have: + +$$\frac{E[\Phi] - E[0]}{L} \approx \frac{1}{2} \Upsilon(L) \bigg (\frac{\Phi}{L}\bigg)^2 + \cdots$$ + +In order to find the ground state under these twisted boundary conditions, we must construct our own variant of the Bose-Hubbard Hamiltonian. Typically you would want to take a peek at the [source code](https://github.com/QuantumKitHub/MPSKitModels.jl/blob/main/src/models/hamiltonians.jl#L379) of `MPSKitModels.jl` to see how these models are defined and tweak it as per your needs. Here we see that applying twisted boundary conditions is equivalent to adding a prefactor of $e^{\pm i\phi}$ in front of the hopping amplitudes. +""" + +function bose_hubbard_model_twisted_bc(elt::Type{<:Number}=ComplexF64, symmetry::Type{<:Sector}=Trivial, + lattice::AbstractLattice=InfiniteChain(1); + cutoff::Integer=5, t=1.0, U=1.0, mu=0.0, phi=0) + + a_pm = a_plusmin(elt, symmetry; cutoff=cutoff) + a_mp = a_minplus(elt, symmetry; cutoff=cutoff) + N = a_number(elt, symmetry; cutoff=cutoff) + + interaction_term = N * (N - id(domain(N))) + + H = @mpoham begin + sum(nearest_neighbours(lattice)) do (i, j) + return -t * (exp(1im * phi) * a_pm{i,j} + exp(1im * -phi) * a_mp{i,j}) + end + + sum(vertices(lattice)) do i + return U / 2 * interaction_term{i} - mu * N{i} + end + end +end + +function superfluid_stiffness_profile(t, mu, D, cutoff, ϵ=1e-4, npoints=11) + phis = range(-ϵ, ϵ, npoints) + energies = zeros(length(phis)) + + Threads.@threads for idx in eachindex(phis) + hamiltonian_twisted = bose_hubbard_model_twisted_bc(cutoff=cutoff, t=t, mu=mu, U=1, phi=phis[idx]) + state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) + state_twisted, _, _ = find_groundstate(state, hamiltonian_twisted, VUMPS(; tol=1e-8, verbosity=0)) + energies[idx] = real(expectation_value(state_twisted, hamiltonian_twisted)) + end + + plot(phis, energies, lw=2, xlabel="Phase twist per site" * L"(\phi)", ylabel="Ground state energy", title="t = $t | μ = $mu | D = $D | cutoff = $cutoff") +end + +superfluid_stiffness_profile(0.2, 0.3, 5, 4) # superfluid + +superfluid_stiffness_profile(0.01, 0.3, 5, 4) # mott insulator + +md""" +Now that we know what phases to expect, we can plot the phase diagram by scanning over a range of parameters. In general, one could do better by performing a bisection algorithm for each chemical potential to determine the value of the hopping parameter at the transition point, however the 1D Bose-Hubbard model may have two transition points at the same chemical potential which makes this a bit cumbersome to implement robustly. Furthermore, we stick to using the quasi-condensate density as an order parameter since extracting the superfluid density accurately requires a more robust scheme to compute second derivatives which takes us away from the focus of this tutorial. +""" + +cutoff, D = 4, 10 +mus = range(0, 0.75, 40) +ts = range(0, 0.3, 40) + +a_op = a_min(cutoff=cutoff) +order_parameters = zeros(length(ts), length(mus)) + +Threads.@threads for (i, j) in collect(Iterators.product(eachindex(mus), eachindex(ts))) + hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff=cutoff, U=1, mu=mus[i], t=ts[j]) + init_state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) + state, _, _ = find_groundstate(init_state, hamiltonian, VUMPS(; tol=1e-8, verbosity=0)) + order_parameters[i, j] = abs(expectation_value(state, 0 => a_op)) +end + +heatmap(ts, mus, order_parameters, xlabel=L"t/U", ylabel=L"\mu/U", title=L"\langle \hat{a}_i \rangle") + +md""" +Although the bond dimension here is quite low, we already see the deformation of the Mott insulator lobes to give way to the well known BKT transition that happens at commensurate density. One can go further and estimate the critical exponents using finite-entanglement scaling procedures on the correlation functions, but these may now be performed with ease using what we have learnt in this tutorial. +""" \ No newline at end of file From cd43827c63a9726cccb0c44dfa5f99fe7ed77a98 Mon Sep 17 00:00:00 2001 From: Akshay Shankar Date: Fri, 31 Oct 2025 15:25:45 +0100 Subject: [PATCH 2/8] formatting --- examples/quantum1d/8.bose-hubbard/main.jl | 122 ++++++++++++---------- 1 file changed, 67 insertions(+), 55 deletions(-) diff --git a/examples/quantum1d/8.bose-hubbard/main.jl b/examples/quantum1d/8.bose-hubbard/main.jl index dd7e2a3c8..665ac9f2d 100644 --- a/examples/quantum1d/8.bose-hubbard/main.jl +++ b/examples/quantum1d/8.bose-hubbard/main.jl @@ -3,7 +3,7 @@ using MPSKit, MPSKitModels, TensorKit using Plots, LaTeXStrings theme(:wong) -default(fontfamily="Computer Modern", label=nothing, dpi=100, framestyle=:box) +default(fontfamily = "Computer Modern", label = nothing, dpi = 100, framestyle = :box) md""" # 1D Bose-Hubbard model @@ -62,21 +62,21 @@ md""" This simply initializes a tensor filled with random entries (check out the documentation for other useful constructors). Next, we need the creation and annihilation operators. While we could construct them from scratch, here we will use `MPSKitModels.jl` instead that has predefined operators and models for most well-known lattice models. """ -a_op = a_min(cutoff=cutoff) # creates a bosonic annihilation operator without any symmetries +a_op = a_min(cutoff = cutoff) # creates a bosonic annihilation operator without any symmetries display(a_op[]) -display((a_op'*a_op)[]) +display((a_op' * a_op)[]) md""" The [] accessor lets us see the underlying array and indeed the operators are exactly what we require. Similarly, the Bose Hubbard model is also predefined (although we will construct our own variant later on). """ -hamiltonian = bose_hubbard_model(InfiniteChain(1), cutoff=cutoff, U=1, mu=0.5, t=0.2) # It is not strictly required to pass InfiniteChain() and is only included for clarity; one may instead pass FiniteChain(N) as well +hamiltonian = bose_hubbard_model(InfiniteChain(1), cutoff = cutoff, U = 1, mu = 0.5, t = 0.2) # It is not strictly required to pass InfiniteChain() and is only included for clarity; one may instead pass FiniteChain(N) as well md""" This has created the Hamiltonian operator as a matrix product operator (MPO) which is a convenient form to use in conjunction with MPS. Finally, the ground state optimization may be performed with either `iDMRG` or `VUMPS`. Both should take similar arguments but it is known that VUMPS is typically more efficient for these systems so we proceed with that. """ -ground_state, _, _ = find_groundstate(initial_state, hamiltonian, VUMPS(tol=1e-6, verbosity=2, maxiter=200)) +ground_state, _, _ = find_groundstate(initial_state, hamiltonian, VUMPS(tol = 1.0e-6, verbosity = 2, maxiter = 200)) println("Energy: ", expectation_value(ground_state, hamiltonian)) md""" @@ -84,30 +84,30 @@ This automatically runs the algorithm until a certain error [measure](https://gi """ function get_ground_state(mu, t, cutoff, D; kwargs...) - hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff=cutoff, U=1, mu=mu, t=t) + hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff = cutoff, U = 1, mu = mu, t = t) state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) state, _, _ = find_groundstate(state, hamiltonian, VUMPS(; kwargs...)) return state end -ground_state = get_ground_state(0.5, 0.01, cutoff, D; tol=1e-6, verbosity=2, maxiter=500) +ground_state = get_ground_state(0.5, 0.01, cutoff, D; tol = 1.0e-6, verbosity = 2, maxiter = 500) md""" Now that we have the state, we may compute observables using the `expectation_value` function. It typically expects a `Pair`, `(i1, i2, .., ik) => op` where `op` is a `TensorMap` or `InfiniteMPO` acting over `k` sites. In case of the Hamiltonian, it is not necessary to specify the indices as it spans the whole lattice. We can now plot the correlation function $\langle \hat{a}^{\dagger}_i \hat{a}_j\rangle$. """ -plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:50), lw=2, xlabel="Site index", ylabel="Correlation function", yscale=:log10) -hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls=:dash, c=:black) +plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:50), lw = 2, xlabel = "Site index", ylabel = "Correlation function", yscale = :log10) +hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls = :dash, c = :black) md""" We see that the correlations drop off exponentially, indicating the existence of a gapped Mott insulating phase. Let us now shift our parameters to probe other phases. """ -ground_state = get_ground_state(0.5, 0.2, cutoff, D; tol=1e-6, verbosity=2, maxiter=500) +ground_state = get_ground_state(0.5, 0.2, cutoff, D; tol = 1.0e-6, verbosity = 2, maxiter = 500) -plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:100), lw=2, xlabel="Site index", ylabel="Correlation function", yscale=:log10, xscale=:log10) -hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls=:dash, c=:black) +plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:100), lw = 2, xlabel = "Site index", ylabel = "Correlation function", yscale = :log10, xscale = :log10) +hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls = :dash, c = :black) md""" In this case, the correlation function drops off algebraically and eventually saturates as $\lim_{i \to \infty}\langle\hat{a}_i^{\dagger} \hat{a}_j\rangle ≈ \langle \hat{a}_i^{\dagger}\rangle \langle \hat{a}_j \rangle = |\langle a_i \rangle|^2 \neq 0$. This is a signature of long-range order and suggests the existence of a Bose-Einstein condensate. However, this is a bit odd since at zero temperature, the Bose Hubbard model is not expected to break any continuous symmetries ($U(1)$ in this case, corresponding to particle number conservation) due to the Mermin-Wagner theorem. The source of this contradiction lies in the fact that the true 1D superfluid ground state is an extended critical phase exhibiting algebraic decay, however, a finite bond-dimension MPS can only capture exponentially decaying correlations. As a result, the finite bond dimension effectively introduces a length scale into the system in a similar manner as finite-size effects. We can see this clearly by increasing the bond dimension. We also see that the correlation length seems to depend algebraically on the bond dimension as expected from finite-entanglement scaling arguments. @@ -119,36 +119,42 @@ mu, t = 0.5, 0.2 states = Vector{InfiniteMPS}(undef, length(Ds)) Threads.@threads for idx in eachindex(Ds) - states[idx] = get_ground_state(mu, t, cutoff, Ds[idx], tol=1e-7, verbosity=1, maxiter=500) + states[idx] = get_ground_state(mu, t, cutoff, Ds[idx], tol = 1.0e-7, verbosity = 1, maxiter = 500) end npoints = 400 two_point_correlation = zeros(length(Ds), npoints) -a_op = a_min(cutoff=cutoff) +a_op = a_min(cutoff = cutoff) Threads.@threads for idx in eachindex(Ds) two_point_correlation[idx, :] .= real.(expectation_value(states[idx], (1, i) => a_op' ⊗ a_op) for i in 1:npoints) end -p = plot(framestyle=:box, ylabel="Correlation function " * L"\langle a_i^{\dagger}a_j \rangle", - xlabel="Distance " * L"|i-j|", xscale=:log10, yscale=:log10, - xticks=([10, 100], ["10", "100"]), - yticks=([0.25, 0.5, 1.0], ["0.25", "0.5", "1.0"])) - -plot!(p, 2:npoints, two_point_correlation[:, 2:end]', - lab="D = " .* string.(permutedims(Ds)), lw=2) - -scatter!(p, Ds, correlation_length.(states), - ylabel="Correlation length", xlabel="Bond dimension", - xscale=:log10, yscale=:log10, - inset=bbox(0.2, 0.51, 0.25, 0.25), - subplot=2, - xticks=(20:10:50, string.(20:10:50)), - yticks=([50, 100], string.([50, 100])), - xlabelfontsize=8, - ylabelfontsize=8, - ylims=[20, 130], - xlims=[15, 60]) +p = plot( + framestyle = :box, ylabel = "Correlation function " * L"\langle a_i^{\dagger}a_j \rangle", + xlabel = "Distance " * L"|i-j|", xscale = :log10, yscale = :log10, + xticks = ([10, 100], ["10", "100"]), + yticks = ([0.25, 0.5, 1.0], ["0.25", "0.5", "1.0"]) +) + +plot!( + p, 2:npoints, two_point_correlation[:, 2:end]', + lab = "D = " .* string.(permutedims(Ds)), lw = 2 +) + +scatter!( + p, Ds, correlation_length.(states), + ylabel = "Correlation length", xlabel = "Bond dimension", + xscale = :log10, yscale = :log10, + inset = bbox(0.2, 0.51, 0.25, 0.25), + subplot = 2, + xticks = (20:10:50, string.(20:10:50)), + yticks = ([50, 100], string.([50, 100])), + xlabelfontsize = 8, + ylabelfontsize = 8, + ylims = [20, 130], + xlims = [15, 60] +) md""" This shows that any finite bond dimension MPS necessarily breaks the symmetry of the system, forming a Bose-Einstein condensate which introduces erraneous long-distance behaviour of correlation functions. In case of finite bond dimension, it is thus reasonable to associate the finite expectation value of the field operator to the 'quasicondensate' density of the system which vanishes as $D \to \infty$. @@ -186,11 +192,15 @@ However, we know that a finite bond dimension MPS introduces a non-zero quasi-co ks = range(-0.05, 0.15, 500) momentum_distribution = map( ((corr, qc),) -> sum( - 2 .* cos.(ks' .* (2:npoints)) .* (corr[2:end] .- qc), dims=1) .+ (corr[1] .- qc), - zip(eachrow(two_point_correlation), - quasicondensate_density)) + 2 .* cos.(ks' .* (2:npoints)) .* (corr[2:end] .- qc), dims = 1 + ) .+ (corr[1] .- qc), + zip( + eachrow(two_point_correlation), + quasicondensate_density + ) +) momentum_distribution = vcat(momentum_distribution...)' -plot(ks, momentum_distribution, lab="D = " .* string.(permutedims(Ds)), lw=1.5, xlabel="Momentum k", ylabel=L"\langle n_k \rangle", ylim=[0, 50]) +plot(ks, momentum_distribution, lab = "D = " .* string.(permutedims(Ds)), lw = 1.5, xlabel = "Momentum k", ylabel = L"\langle n_k \rangle", ylim = [0, 50]) md""" We see that the density seems to peak around $k=0$, this time seemingly becoming more prominent as $D \to \infty$ which seems to suggest again that there is a condensate. However, going by the Penrose-Onsager criterion, the existence of a condensate can be quantified by requiring the leading eigenvalue of the single particle density matrix (i.e, $\langle \hat{n}_{k=0}\rangle = \sum_j \langle \hat{a}_j^{\dagger} \hat{a}_0\rangle$) to diverge as $O(N)$ in the thermodynamic limit. In this case, since the correlations decay as a power law, there is naturally a divergence at low momenta. But this does not imply the existence of a condensate since the order of divergence is much weaker. However, this does indicate the remnants of some kind of condensation in the 1D model despite the quantum fluctuations, leading to the practical utility of defining the concept of a quasicondensate where there is still a notion of phase coherence over short distances. @@ -202,38 +212,40 @@ $$\frac{E[\Phi] - E[0]}{L} \approx \frac{1}{2} \Upsilon(L) \bigg (\frac{\Phi}{L} In order to find the ground state under these twisted boundary conditions, we must construct our own variant of the Bose-Hubbard Hamiltonian. Typically you would want to take a peek at the [source code](https://github.com/QuantumKitHub/MPSKitModels.jl/blob/main/src/models/hamiltonians.jl#L379) of `MPSKitModels.jl` to see how these models are defined and tweak it as per your needs. Here we see that applying twisted boundary conditions is equivalent to adding a prefactor of $e^{\pm i\phi}$ in front of the hopping amplitudes. """ -function bose_hubbard_model_twisted_bc(elt::Type{<:Number}=ComplexF64, symmetry::Type{<:Sector}=Trivial, - lattice::AbstractLattice=InfiniteChain(1); - cutoff::Integer=5, t=1.0, U=1.0, mu=0.0, phi=0) +function bose_hubbard_model_twisted_bc( + elt::Type{<:Number} = ComplexF64, symmetry::Type{<:Sector} = Trivial, + lattice::AbstractLattice = InfiniteChain(1); + cutoff::Integer = 5, t = 1.0, U = 1.0, mu = 0.0, phi = 0 + ) - a_pm = a_plusmin(elt, symmetry; cutoff=cutoff) - a_mp = a_minplus(elt, symmetry; cutoff=cutoff) - N = a_number(elt, symmetry; cutoff=cutoff) + a_pm = a_plusmin(elt, symmetry; cutoff = cutoff) + a_mp = a_minplus(elt, symmetry; cutoff = cutoff) + N = a_number(elt, symmetry; cutoff = cutoff) interaction_term = N * (N - id(domain(N))) - H = @mpoham begin + return H = @mpoham begin sum(nearest_neighbours(lattice)) do (i, j) - return -t * (exp(1im * phi) * a_pm{i,j} + exp(1im * -phi) * a_mp{i,j}) + return -t * (exp(1im * phi) * a_pm{i, j} + exp(1im * -phi) * a_mp{i, j}) end + - sum(vertices(lattice)) do i + sum(vertices(lattice)) do i return U / 2 * interaction_term{i} - mu * N{i} end end end -function superfluid_stiffness_profile(t, mu, D, cutoff, ϵ=1e-4, npoints=11) +function superfluid_stiffness_profile(t, mu, D, cutoff, ϵ = 1.0e-4, npoints = 11) phis = range(-ϵ, ϵ, npoints) energies = zeros(length(phis)) Threads.@threads for idx in eachindex(phis) - hamiltonian_twisted = bose_hubbard_model_twisted_bc(cutoff=cutoff, t=t, mu=mu, U=1, phi=phis[idx]) + hamiltonian_twisted = bose_hubbard_model_twisted_bc(cutoff = cutoff, t = t, mu = mu, U = 1, phi = phis[idx]) state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) - state_twisted, _, _ = find_groundstate(state, hamiltonian_twisted, VUMPS(; tol=1e-8, verbosity=0)) + state_twisted, _, _ = find_groundstate(state, hamiltonian_twisted, VUMPS(; tol = 1.0e-8, verbosity = 0)) energies[idx] = real(expectation_value(state_twisted, hamiltonian_twisted)) end - plot(phis, energies, lw=2, xlabel="Phase twist per site" * L"(\phi)", ylabel="Ground state energy", title="t = $t | μ = $mu | D = $D | cutoff = $cutoff") + return plot(phis, energies, lw = 2, xlabel = "Phase twist per site" * L"(\phi)", ylabel = "Ground state energy", title = "t = $t | μ = $mu | D = $D | cutoff = $cutoff") end superfluid_stiffness_profile(0.2, 0.3, 5, 4) # superfluid @@ -248,18 +260,18 @@ cutoff, D = 4, 10 mus = range(0, 0.75, 40) ts = range(0, 0.3, 40) -a_op = a_min(cutoff=cutoff) +a_op = a_min(cutoff = cutoff) order_parameters = zeros(length(ts), length(mus)) Threads.@threads for (i, j) in collect(Iterators.product(eachindex(mus), eachindex(ts))) - hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff=cutoff, U=1, mu=mus[i], t=ts[j]) + hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff = cutoff, U = 1, mu = mus[i], t = ts[j]) init_state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) - state, _, _ = find_groundstate(init_state, hamiltonian, VUMPS(; tol=1e-8, verbosity=0)) + state, _, _ = find_groundstate(init_state, hamiltonian, VUMPS(; tol = 1.0e-8, verbosity = 0)) order_parameters[i, j] = abs(expectation_value(state, 0 => a_op)) end -heatmap(ts, mus, order_parameters, xlabel=L"t/U", ylabel=L"\mu/U", title=L"\langle \hat{a}_i \rangle") +heatmap(ts, mus, order_parameters, xlabel = L"t/U", ylabel = L"\mu/U", title = L"\langle \hat{a}_i \rangle") md""" Although the bond dimension here is quite low, we already see the deformation of the Mott insulator lobes to give way to the well known BKT transition that happens at commensurate density. One can go further and estimate the critical exponents using finite-entanglement scaling procedures on the correlation functions, but these may now be performed with ease using what we have learnt in this tutorial. -""" \ No newline at end of file +""" From ef8a2e61279740908dca118c13c3584b255a7d45 Mon Sep 17 00:00:00 2001 From: Akshay Shankar Date: Fri, 31 Oct 2025 15:43:05 +0100 Subject: [PATCH 3/8] update notebooks after formatting --- .../quantum1d/8.bose-hubbard/index.md | 158 ++++++++++-------- .../quantum1d/8.bose-hubbard/main.ipynb | 120 +++++++------ 2 files changed, 151 insertions(+), 127 deletions(-) diff --git a/docs/src/examples/quantum1d/8.bose-hubbard/index.md b/docs/src/examples/quantum1d/8.bose-hubbard/index.md index 1e1a1af2e..3396719f3 100644 --- a/docs/src/examples/quantum1d/8.bose-hubbard/index.md +++ b/docs/src/examples/quantum1d/8.bose-hubbard/index.md @@ -12,7 +12,7 @@ using MPSKit, MPSKitModels, TensorKit using Plots, LaTeXStrings theme(:wong) -default(fontfamily="Computer Modern", label=nothing, dpi=100, framestyle=:box) +default(fontfamily = "Computer Modern", label = nothing, dpi = 100, framestyle = :box) ```` # 1D Bose-Hubbard model @@ -71,8 +71,8 @@ initial_state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) ```` single site InfiniteMPS: │ ⋮ -│ C[1]: TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 1, Vector{ComplexF64}}(ComplexF64[0.97751+0.0im, 0.174631+0.0125955im, 0.104701+0.000310339im, 0.0116397-0.00627461im, 0.0464636+0.000250548im, 0.0+0.0im, 0.0143524+0.0im, 0.00360949+0.00128115im, -0.00368336-0.0047228im, 0.000692314+0.00211969im, 0.0+0.0im, 0.0+0.0im, 0.00655329+0.0im, 0.00552591+0.00129082im, -0.00303456+0.00143664im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0102319+0.0im, -0.00205334+0.00505879im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.00640026+0.0im], ℂ^5 ← ℂ^5) -├── AL[1]: TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[0.390809+0.214646im, 0.0489818+0.0751624im, 0.0497281+0.0204177im, 0.0263021-0.000245521im, 0.0120706+0.0101434im, 0.211349+0.301151im, 0.0250289+0.0560838im, 0.0305354+0.027863im, 0.0167841+3.42096e-5im, 0.0149195+0.010686im, 0.324269+0.45im, 0.0444069+0.0705503im, 0.0396325+0.0521653im, 0.00128175+0.00922583im, 0.0179868+0.0194204im, 0.29577+0.274281im, 0.0536787+0.0464249im, 0.0358399+0.0188768im, 0.0105704-0.00635067im, 0.00311747+0.0202148im, 0.234359+0.307432im, 0.0459379+0.0535219im, 0.0237564+0.0496054im, -0.00216268+0.0275897im, 0.00395079+0.000283255im, -0.24849+0.340759im, -0.0387098-0.00943688im, -0.066458+0.00845679im, -0.0827909-0.00877386im, 0.0442285+0.0196632im, 0.455119-0.222154im, 0.187925-0.0124478im, 0.0290096-0.0172291im, -0.0627177-0.0544896im, 0.0249878+0.0111796im, -0.258925-0.366019im, -0.0495616+0.00700246im, -0.0334674-0.0427697im, -0.0235555+0.00522338im, 0.0129344-0.023551im, 0.370142+0.3491im, 0.0503108+0.0969619im, 0.0315893+0.0386184im, 0.00266811+0.00258244im, 0.0500739+0.0254822im, 0.0565483+0.077467im, 0.0203963+0.0620186im, -0.00282502-0.0492105im, 0.0218008-0.057548im, -0.00380262+0.0459115im, 0.216276+0.596959im, 0.103859+0.0332121im, 0.031629+0.105885im, 0.0856072+0.0532032im, -0.0631615+0.0551744im, 0.0300362-0.137117im, -0.086399-0.0116169im, 0.0315047+0.0466518im, 0.01568+0.0662402im, 0.000724708-0.0406109im, 0.234667-0.308019im, 0.167079-0.192054im, 0.00466488-0.0928214im, 0.00539646-0.0186621im, 0.00962596-0.0344942im, -0.456697-0.139859im, -0.0242026-0.0025607im, -0.038756-0.014274im, 0.0144336+0.0956192im, -0.0437668-0.0752814im, 0.107807+0.107598im, -0.00732072+0.0396839im, 0.102896+0.0301685im, 0.0688301-0.0658399im, 0.0317364-0.00169938im, 0.0520469-0.259208im, -0.0400713-0.0759im, -0.0490438+0.0512775im, -0.141529+0.0757342im, 0.0237932-0.0950623im, -0.109863-0.118982im, 0.0226764-0.193153im, -0.094942-0.070543im, -0.0890496-0.0436445im, 0.0274223+0.0547519im, -0.039393-0.218089im, 0.0163183+0.21468im, 0.0132428-0.0222767im, 0.0425839-0.114359im, 0.0249832+0.0736523im, -0.121223-0.0674908im, -0.0640196+0.0505719im, -0.0662815+0.026577im, -0.0533572-0.126597im, 0.0568246+0.0547217im, 0.169584+0.743324im, -0.00516459+0.0778608im, -0.0791716+0.0630091im, -0.0966301+0.00510259im, 0.0637203+0.0805271im, 0.140845+0.171094im, 0.021739-0.177908im, -0.0762735+0.11312im, -0.191275+0.0465106im, 0.0304639-0.0648004im, 0.229036+0.518716im, 0.0905234+0.039573im, -0.0974247-0.0245568im, 0.0249835-0.092148im, -0.0546524+0.101555im, -0.222658-0.0990489im, -0.198631+0.169606im, -0.0657467+0.0130075im, 0.0849531-0.0726054im, -0.0655968+0.0338763im, -0.0709195-0.353206im, -0.119811-0.0898848im, -0.0932307+0.106653im, -0.10785-0.0538011im, 0.0929571+0.0313938im, 0.132025-0.162429im, 0.0960656-0.142876im, -0.0589741-0.111375im, -0.14054-0.142388im, 0.0681092+0.109553im], (ℂ^5 ⊗ ℂ^5) ← ℂ^5) +│ C[1]: TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 1, Vector{ComplexF64}}(ComplexF64[0.980773+0.0im, 0.146227-0.0278782im, 0.105333-0.010143im, 0.0620678+0.00292624im, 0.0241261-0.00222812im, 0.0+0.0im, 0.0101049+0.0im, 0.00743715-0.00597138im, -0.0011665+0.00106208im, -0.00256436-0.000271479im, 0.0+0.0im, 0.0+0.0im, 0.00580994+0.0im, -0.00381868-0.00101898im, -0.00128419+0.00100647im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.00397274+0.0im, 0.00234497-0.000379053im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.00198522+0.0im], ℂ^5 ← ℂ^5) +├── AL[1]: TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[0.375885+0.305768im, 0.0768109+0.0330615im, 0.0460733+0.0360179im, 0.0240341+0.025383im, 0.00804042+0.0057914im, 0.213348+0.271799im, 0.0517906+0.0297947im, 0.0138421+0.0188126im, 0.00857498+0.0239495im, 0.00590194+0.0070995im, 0.393411+0.386699im, 0.0781762+0.0467087im, 0.0547209+0.0259146im, 0.0239884+0.026463im, 0.0118271+0.00963569im, 0.242906+0.273226im, 0.0539362+0.037471im, 0.0475866+0.0266459im, 0.0156312+0.0217225im, 0.0105492+0.00484183im, 0.195707+0.347509im, 0.0485176+0.057361im, 0.0399442+0.0476088im, 0.0108898+0.0292134im, 0.00742594+0.00871371im, -0.405789-0.094101im, -0.00619161+0.0455026im, -0.0225547-0.107661im, -0.0576217-0.0382695im, -0.0256341-0.0202461im, 0.480858+0.520426im, 0.000291008+0.091172im, 0.139655+0.119484im, 0.0651289-0.0131326im, 0.017084+0.0264884im, -0.21137-0.149131im, -0.0469705-0.104014im, -0.154879-0.00320825im, -0.0192648-0.007215im, -0.0348013-0.0153511im, 0.0504188+0.0087985im, -0.0921246-0.012542im, -0.120437+0.0817629im, 0.0133867-0.0348062im, -0.0271514-0.0101404im, 0.294447+0.113036im, 0.0607909-0.0783091im, -0.0267377-0.0718419im, 0.0296402-0.0108425im, -0.00573401-0.0122223im, 0.0103164-0.142274im, -0.0646881-0.0987699im, -0.0188412+0.0287894im, 0.0265472-0.0232285im, -0.00197098+0.00516289im, 0.0692401-0.176396im, 0.0426622-0.0273936im, -0.00973839-0.0308723im, 0.0120297+0.0305825im, 0.0208432-0.00767607im, -0.251962-0.320331im, -0.0961427-0.00484129im, 0.0287567-0.0230259im, 0.0149031-0.0233169im, 0.0280275-0.00709338im, 0.371239+0.719013im, 0.114536+0.131076im, 0.123631+0.02352im, -0.0317699+0.0460085im, 0.00597557+0.0348219im, -0.125585+0.0936419im, -0.0818155+0.0358405im, -0.0724073+0.022764im, -0.0116091-0.00616551im, 0.00110864+0.00807954im, 0.40961+0.365138im, -0.128697+0.130384im, 0.118639+0.107234im, 0.0506963-0.0355463im, 0.0413307+0.0240077im, -0.125027+0.164525im, -0.113642-0.0873257im, -0.137262-0.00708175im, 0.0456722+0.0461297im, 0.00611336-0.0213785im, -0.39764-0.209872im, -0.00809221+0.0215245im, -0.026733-0.00806058im, -0.0233946-0.0571397im, 0.0202147-0.0273948im, -0.0566267+0.0987188im, -0.0429491-0.0913588im, -0.179772-0.115045im, 0.0264671+0.0644071im, 0.0164148+0.00709543im, 0.353547-0.317267im, 0.108454-0.00632667im, 0.114669-0.10561im, 0.0304626+0.0034387im, 0.041428-0.00870923im, -0.113545-0.159728im, 0.0926693-0.0340042im, 0.006828-0.0369635im, -0.121696+0.0542231im, -0.0053801+0.0419627im, -0.398247-0.168103im, -0.14182+0.206577im, 0.0674782+0.207936im, -0.0729256-0.175661im, -0.0606546+0.0189888im, 0.0538273+0.214954im, -0.183298-0.00539617im, -0.1011+0.0842453im, -0.052954+0.0404571im, -0.0774749+0.0641526im, -0.0490675+0.214829im, -0.00844192+0.151745im, 0.105707+0.130237im, 0.0499562-0.0486145im, 0.0121478-0.0154518im, 0.572407-0.0611613im, -0.0212294-0.128364im, -0.00861733-0.109909im, -0.00619022-0.07167im, -0.0857339-0.0245475im], (ℂ^5 ⊗ ℂ^5) ← ℂ^5) │ ⋮ ```` @@ -80,15 +80,15 @@ single site InfiniteMPS: This simply initializes a tensor filled with random entries (check out the documentation for other useful constructors). Next, we need the creation and annihilation operators. While we could construct them from scratch, here we will use `MPSKitModels.jl` instead that has predefined operators and models for most well-known lattice models. ````julia -a_op = a_min(cutoff=cutoff) # creates a bosonic annihilation operator without any symmetries +a_op = a_min(cutoff = cutoff) # creates a bosonic annihilation operator without any symmetries display(a_op[]) -display((a_op'*a_op)[]) +display((a_op' * a_op)[]) ```` The [] accessor lets us see the underlying array and indeed the operators are exactly what we require. Similarly, the Bose Hubbard model is also predefined (although we will construct our own variant later on). ````julia -hamiltonian = bose_hubbard_model(InfiniteChain(1), cutoff=cutoff, U=1, mu=0.5, t=0.2) # It is not strictly required to pass InfiniteChain() and is only included for clarity; one may instead pass FiniteChain(N) as well +hamiltonian = bose_hubbard_model(InfiniteChain(1), cutoff = cutoff, U = 1, mu = 0.5, t = 0.2) # It is not strictly required to pass InfiniteChain() and is only included for clarity; one may instead pass FiniteChain(N) as well ```` ```` @@ -102,14 +102,14 @@ single site InfiniteMPOHamiltonian{MPSKit.JordanMPOTensor{ComplexF64, TensorKit. This has created the Hamiltonian operator as a matrix product operator (MPO) which is a convenient form to use in conjunction with MPS. Finally, the ground state optimization may be performed with either `iDMRG` or `VUMPS`. Both should take similar arguments but it is known that VUMPS is typically more efficient for these systems so we proceed with that. ````julia -ground_state, _, _ = find_groundstate(initial_state, hamiltonian, VUMPS(tol=1e-6, verbosity=2, maxiter=200)) +ground_state, _, _ = find_groundstate(initial_state, hamiltonian, VUMPS(tol = 1.0e-6, verbosity = 2, maxiter = 200)) println("Energy: ", expectation_value(ground_state, hamiltonian)) ```` ```` -[ Info: VUMPS init: obj = +3.785371201331e-01 err = 5.9503e-01 -[ Info: VUMPS conv 48: obj = -6.757777651160e-01 err = 8.4647387076e-07 time = 13.78 sec -Energy: -0.6757777651160398 - 5.389840049968564e-17im +[ Info: VUMPS init: obj = +3.877905708364e-01 err = 6.1967e-01 +[ Info: VUMPS conv 33: obj = -6.757777651162e-01 err = 8.3203601611e-07 time = 9.16 sec +Energy: -0.675777765116157 + 6.114900252818245e-17im ```` @@ -117,21 +117,21 @@ This automatically runs the algorithm until a certain error [measure](https://gi ````julia function get_ground_state(mu, t, cutoff, D; kwargs...) - hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff=cutoff, U=1, mu=mu, t=t) + hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff = cutoff, U = 1, mu = mu, t = t) state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) state, _, _ = find_groundstate(state, hamiltonian, VUMPS(; kwargs...)) return state end -ground_state = get_ground_state(0.5, 0.01, cutoff, D; tol=1e-6, verbosity=2, maxiter=500) +ground_state = get_ground_state(0.5, 0.01, cutoff, D; tol = 1.0e-6, verbosity = 2, maxiter = 500) ```` ```` single site InfiniteMPS: │ ⋮ -│ C[1]: TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 1, Vector{ComplexF64}}(ComplexF64[0.54949+0.0im, 0.369229+0.0647129im, 0.178404+0.141066im, -0.315274+0.540671im, -0.123933+0.313108im, 0.0+0.0im, 0.0146972+0.0im, 0.00310867+0.00509383im, -0.00170132+0.0120718im, -0.00584434+0.00729265im, 0.0+0.0im, 0.0+0.0im, 0.0120375+0.0im, -0.00309885-0.000910003im, -0.00472919-0.00552589im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 3.69614e-6+0.0im, 1.57727e-6+3.97113e-7im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 3.35037e-6+0.0im], ℂ^5 ← ℂ^5) -├── AL[1]: TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[0.334572+0.0477944im, 0.0555237+0.176413im, 0.166278+0.112267im, -0.33895+0.132832im, -0.124122-0.025603im, 0.243524+0.200797im, 0.12627+0.147902im, 0.0280863+0.120709im, -0.321951+0.11476im, -0.154783+0.0931863im, 0.188335+0.150819im, 0.0985741+0.105361im, 0.0296448+0.0819545im, -0.241768+0.094309im, -0.120223+0.0771641im, 0.124161+0.204734im, 0.0575478+0.149159im, -0.0126871+0.0996538im, -0.269509+0.00302196im, -0.140838+0.0244051im, 0.000861337+0.00207717im, -6.42551e-5-0.00173344im, -0.00105528+0.00179777im, 0.000764672-0.000602481im, 0.00189714+0.00127471im, -0.30991+0.121537im, -0.0319424-0.144377im, -0.231995-0.0529414im, 0.241863-0.170549im, 0.0820029+0.0758833im, 0.168514+0.0935542im, 0.117687+0.0982952im, 0.0270232+0.081966im, -0.202623+0.122908im, -0.104325+0.0771085im, -0.156916-0.161661im, -0.101529-0.147743im, -0.000801824-0.117147im, 0.265459-0.0682365im, 0.135124-0.0486043im, 0.200206+0.268495im, 0.104333+0.207532im, -0.00332133+0.137411im, -0.382703+0.0442051im, -0.202173+0.0530631im, 0.00151036+0.00282271im, -0.000333829-0.00244654im, -0.00124394+0.002695im, 0.00099978-0.000960587im, 0.00286522+0.00154189im, 0.0811625-0.136874im, 0.286116+0.0340031im, 0.039079-0.115171im, -0.0396969+0.373233im, -0.160445+0.225487im, 0.12465+0.00265451im, 0.0853333+0.0120395im, 0.0583849+0.0470931im, -0.0738299+0.118299im, -0.0280568+0.0600721im, -0.386529+0.140124im, -0.270394+0.0470919im, -0.154767-0.0447004im, 0.0837792-0.458016im, 0.0072124-0.2535im, -0.133538+0.0105836im, -0.084408-0.00812938im, -0.0475752-0.0328252im, 0.0654166-0.130799im, 0.0203034-0.0724504im, -0.00122441-9.18638e-5im, 0.000860446+0.000395458im, -0.000616071-0.000956551im, 0.000105274+0.000520327im, -0.00108792+0.000611516im, 0.0298865-0.142912im, -0.217501+0.180601im, 0.175454-0.030206im, -0.13986-0.160109im, -0.0507848-0.328219im, 0.0525566-0.337304im, 0.0799062-0.22303im, 0.104006-0.0939662im, 0.304141+0.249739im, 0.180722+0.109822im, -0.0592788-0.0217604im, -0.0322979-0.0333287im, 0.013151-0.0260831im, 0.0586576-0.0412738im, 0.01969-0.0305178im, -0.134283+0.272541im, -0.121871+0.168365im, -0.113483+0.0536486im, -0.192178-0.28808im, -0.126192-0.138043im, -0.00164022+0.00237254im, 0.00182384-0.00127475im, -0.00267205+0.000112583im, 0.00118167+0.000405385im, -5.59629e-5+0.00293166im, -0.0899756+0.542342im, 0.284935+0.0481453im, -0.359211+0.11573im, -0.188399+0.0301327im, -0.201842+0.371847im, 0.061353-0.172194im, 0.0567285-0.113556im, 0.0525202-0.0403404im, 0.142324+0.156115im, 0.0951062+0.0784393im, 0.122461-0.0716501im, 0.0975112-0.0334494im, 0.0685959+0.0160209im, -0.00270279+0.165052im, 0.0093008+0.0828786im, -0.141691+0.0383132im, -0.101863+0.00900411im, -0.0553673-0.0233559im, 0.0437576-0.163564im, 0.0112594-0.0913605im, -0.00139728+0.000160781im, 0.00105287+0.000259349im, -0.0008995-0.000943779im, 0.000230458+0.000563954im, -0.00109281+0.000920657im], (ℂ^5 ⊗ ℂ^5) ← ℂ^5) +│ C[1]: TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 1, Vector{ComplexF64}}(ComplexF64[0.259222+0.0im, 0.246018-0.345952im, 0.610738+0.473339im, -0.148954-0.282818im, -0.211517+0.0836342im, 0.0+0.0im, 0.0136571+0.0im, -0.00924763+0.0239744im, 0.0102051-0.00219558im, 0.00373234-0.000401849im, 0.0+0.0im, 0.0+0.0im, 0.0227273+0.0im, -0.00604063-0.00171163im, 0.000546388+0.0105763im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 3.50256e-6+0.0im, 2.80657e-6-1.34091e-6im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 4.26623e-6+0.0im], ℂ^5 ← ℂ^5) +├── AL[1]: TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[-0.087641+0.413009im, 0.214853+0.0747518im, 0.405403+0.3661im, 0.104705-0.379309im, 0.00466853-0.282206im, 0.0435006+0.0621785im, 0.122043-0.0029551im, 0.0130668+0.192364im, 0.0341526-0.0779833im, -0.0402589-0.0307335im, 0.0439795+0.0220837im, 0.077068-0.0405486im, 0.0546722+0.118445im, 0.00333692-0.0560385im, -0.0295748-0.00984742im, 0.0631303-0.0715436im, -0.0365518-0.15324im, 0.283061-0.0553091im, -0.115614-0.0282307im, -0.0282239+0.0789666im, 0.00187513-0.00134459im, 0.000525481-0.000776215im, -0.00180619-3.04145e-5im, -0.000788208+0.000901196im, -0.00130384+0.00111009im, 0.237015-0.0571072im, 0.166803-0.0622261im, -0.197895+0.172091im, -0.0499605+0.0113682im, -0.192109+0.0651415im, -0.021684+0.115554im, 0.121816+0.150109im, -0.260216+0.190125im, 0.126714-0.0298417im, -0.0158932-0.0920934im, 0.128831+0.120805im, 0.271616-0.0633097im, 0.102193+0.525133im, 0.0494367-0.218492im, -0.157755-0.0575024im, -0.0885049+0.036656im, -0.0422998+0.154893im, -0.276699-0.0966281im, 0.0859221+0.0802997im, 0.061836-0.0591176im, -0.00232231+0.00025207im, -0.000850451+0.000423277im, 0.00157813+0.000926226im, 0.0011447-0.000402135im, 0.0017028-0.000329999im, -0.19475+0.104128im, -0.127458+0.218784im, -0.197587-0.269976im, 0.0936674+0.0793825im, 0.154231-0.116272im, 0.183613+0.0437624im, 0.241292-0.214837im, 0.38749+0.468974im, -0.058015-0.238825im, -0.164866+0.0298767im, -0.012259+0.0545234im, 0.0492084+0.0688762im, -0.140159+0.0922497im, 0.0632274-0.0161014im, -0.014811-0.0562813im, -0.071732-0.0151331im, -0.089349+0.0791879im, -0.134631-0.168392im, 0.0228751+0.0856916im, 0.0637356-0.0103791im, -0.00155239-0.000769454im, -0.000700543-7.51919e-5im, 0.000616444+0.0012093im, 0.000875375+0.00020615im, 0.00119626+0.000473548im, -0.210309+0.314047im, -0.146014+0.176102im, 0.267523-0.412974im, -0.0061229-0.115727im, 0.153847-0.254992im, -0.0890551+0.0141788im, -0.0472752+0.120377im, -0.215729-0.106406im, 0.0743765+0.0813334im, 0.0811442-0.0391826im, 0.0468722-0.00895011im, 0.0236013-0.0654369im, 0.152916+0.0731364im, -0.0483316-0.0501106im, -0.050022+0.038829im, -0.0458048+0.137925im, 0.13997+0.192221im, -0.360059+0.239392im, 0.176409-0.0288231im, -0.00697949-0.127442im, -0.00182017+0.00294987im, -0.000279671+0.00138125im, 0.00252286-0.00100342im, 0.000572153-0.00170456im, 0.00116508-0.00229521im, -0.123356-0.0834066im, -0.205016+0.0435919im, -0.0418691-0.332833im, -0.0181937+0.139342im, 0.12447+0.0348314im, -0.0212312-0.0617839im, -0.0831529-0.0234744im, 0.0636293-0.159741im, -0.0468686+0.0599996im, 0.0460042+0.0572296im, 0.158391+0.0558601im, 0.241495-0.160082im, 0.252877+0.417002im, -0.0178357-0.198992im, -0.12711-0.000907711im, 0.11528-0.0948115im, -0.0221029-0.24102im, 0.439961-0.0282941im, -0.17246-0.066882im, -0.0626087+0.113866im, 0.00313034-0.00165794im, 0.000960262-0.00107186im, -0.00274142-0.000420511im, -0.0013851+0.00120634im, -0.0022133+0.00141929im], (ℂ^5 ⊗ ℂ^5) ← ℂ^5) │ ⋮ ```` @@ -139,21 +139,21 @@ single site InfiniteMPS: Now that we have the state, we may compute observables using the `expectation_value` function. It typically expects a `Pair`, `(i1, i2, .., ik) => op` where `op` is a `TensorMap` or `InfiniteMPO` acting over `k` sites. In case of the Hamiltonian, it is not necessary to specify the indices as it spans the whole lattice. We can now plot the correlation function $\langle \hat{a}^{\dagger}_i \hat{a}_j\rangle$. ````julia -plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:50), lw=2, xlabel="Site index", ylabel="Correlation function", yscale=:log10) -hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls=:dash, c=:black) +plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:50), lw = 2, xlabel = "Site index", ylabel = "Correlation function", yscale = :log10) +hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls = :dash, c = :black) ```` ```@raw html - + ``` We see that the correlations drop off exponentially, indicating the existence of a gapped Mott insulating phase. Let us now shift our parameters to probe other phases. ````julia -ground_state = get_ground_state(0.5, 0.2, cutoff, D; tol=1e-6, verbosity=2, maxiter=500) +ground_state = get_ground_state(0.5, 0.2, cutoff, D; tol = 1.0e-6, verbosity = 2, maxiter = 500) -plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:100), lw=2, xlabel="Site index", ylabel="Correlation function", yscale=:log10, xscale=:log10) -hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls=:dash, c=:black) +plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:100), lw = 2, xlabel = "Site index", ylabel = "Correlation function", yscale = :log10, xscale = :log10) +hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls = :dash, c = :black) ```` ```@raw html @@ -169,40 +169,46 @@ mu, t = 0.5, 0.2 states = Vector{InfiniteMPS}(undef, length(Ds)) Threads.@threads for idx in eachindex(Ds) - states[idx] = get_ground_state(mu, t, cutoff, Ds[idx], tol=1e-7, verbosity=1, maxiter=500) + states[idx] = get_ground_state(mu, t, cutoff, Ds[idx], tol = 1.0e-7, verbosity = 1, maxiter = 500) end npoints = 400 two_point_correlation = zeros(length(Ds), npoints) -a_op = a_min(cutoff=cutoff) +a_op = a_min(cutoff = cutoff) Threads.@threads for idx in eachindex(Ds) two_point_correlation[idx, :] .= real.(expectation_value(states[idx], (1, i) => a_op' ⊗ a_op) for i in 1:npoints) end -p = plot(framestyle=:box, ylabel="Correlation function " * L"\langle a_i^{\dagger}a_j \rangle", - xlabel="Distance " * L"|i-j|", xscale=:log10, yscale=:log10, - xticks=([10, 100], ["10", "100"]), - yticks=([0.25, 0.5, 1.0], ["0.25", "0.5", "1.0"])) - -plot!(p, 2:npoints, two_point_correlation[:, 2:end]', - lab="D = " .* string.(permutedims(Ds)), lw=2) - -scatter!(p, Ds, correlation_length.(states), - ylabel="Correlation length", xlabel="Bond dimension", - xscale=:log10, yscale=:log10, - inset=bbox(0.2, 0.51, 0.25, 0.25), - subplot=2, - xticks=(20:10:50, string.(20:10:50)), - yticks=([50, 100], string.([50, 100])), - xlabelfontsize=8, - ylabelfontsize=8, - ylims=[20, 130], - xlims=[15, 60]) +p = plot( + framestyle = :box, ylabel = "Correlation function " * L"\langle a_i^{\dagger}a_j \rangle", + xlabel = "Distance " * L"|i-j|", xscale = :log10, yscale = :log10, + xticks = ([10, 100], ["10", "100"]), + yticks = ([0.25, 0.5, 1.0], ["0.25", "0.5", "1.0"]) +) + +plot!( + p, 2:npoints, two_point_correlation[:, 2:end]', + lab = "D = " .* string.(permutedims(Ds)), lw = 2 +) + +scatter!( + p, Ds, correlation_length.(states), + ylabel = "Correlation length", xlabel = "Bond dimension", + xscale = :log10, yscale = :log10, + inset = bbox(0.2, 0.51, 0.25, 0.25), + subplot = 2, + xticks = (20:10:50, string.(20:10:50)), + yticks = ([50, 100], string.([50, 100])), + xlabelfontsize = 8, + ylabelfontsize = 8, + ylims = [20, 130], + xlims = [15, 60] +) ```` ```@raw html - + ``` This shows that any finite bond dimension MPS necessarily breaks the symmetry of the system, forming a Bose-Einstein condensate which introduces erraneous long-distance behaviour of correlation functions. In case of finite bond dimension, it is thus reasonable to associate the finite expectation value of the field operator to the 'quasicondensate' density of the system which vanishes as $D \to \infty$. @@ -213,13 +219,13 @@ quasicondensate_density = map(state -> abs2(expectation_value(state, (0,) => a_o ```` 7-element Vector{Float64}: - 0.3097427784562408 - 0.2881477510455394 - 0.2702002091214058 - 0.25712728499126414 - 0.24685384625099716 - 0.23539796049808664 - 0.22799690146561793 + 0.30974277582170107 + 0.28814776444605805 + 0.27020016270456915 + 0.257127287513971 + 0.24685386110383362 + 0.23539802368529517 + 0.22799655546142042 ```` We may now also visualize the momentum distribution function, which is obtained as the Fourier transform of the single-particle density matrix. Starting from the definition of the momentum occupation operator: @@ -251,15 +257,19 @@ However, we know that a finite bond dimension MPS introduces a non-zero quasi-co ks = range(-0.05, 0.15, 500) momentum_distribution = map( ((corr, qc),) -> sum( - 2 .* cos.(ks' .* (2:npoints)) .* (corr[2:end] .- qc), dims=1) .+ (corr[1] .- qc), - zip(eachrow(two_point_correlation), - quasicondensate_density)) + 2 .* cos.(ks' .* (2:npoints)) .* (corr[2:end] .- qc), dims = 1 + ) .+ (corr[1] .- qc), + zip( + eachrow(two_point_correlation), + quasicondensate_density + ) +) momentum_distribution = vcat(momentum_distribution...)' -plot(ks, momentum_distribution, lab="D = " .* string.(permutedims(Ds)), lw=1.5, xlabel="Momentum k", ylabel=L"\langle n_k \rangle", ylim=[0, 50]) +plot(ks, momentum_distribution, lab = "D = " .* string.(permutedims(Ds)), lw = 1.5, xlabel = "Momentum k", ylabel = L"\langle n_k \rangle", ylim = [0, 50]) ```` ```@raw html - + ``` We see that the density seems to peak around $k=0$, this time seemingly becoming more prominent as $D \to \infty$ which seems to suggest again that there is a condensate. However, going by the Penrose-Onsager criterion, the existence of a condensate can be quantified by requiring the leading eigenvalue of the single particle density matrix (i.e, $\langle \hat{n}_{k=0}\rangle = \sum_j \langle \hat{a}_j^{\dagger} \hat{a}_0\rangle$) to diverge as $O(N)$ in the thermodynamic limit. In this case, since the correlations decay as a power law, there is naturally a divergence at low momenta. But this does not imply the existence of a condensate since the order of divergence is much weaker. However, this does indicate the remnants of some kind of condensation in the 1D model despite the quantum fluctuations, leading to the practical utility of defining the concept of a quasicondensate where there is still a notion of phase coherence over short distances. @@ -271,38 +281,40 @@ $$\frac{E[\Phi] - E[0]}{L} \approx \frac{1}{2} \Upsilon(L) \bigg (\frac{\Phi}{L} In order to find the ground state under these twisted boundary conditions, we must construct our own variant of the Bose-Hubbard Hamiltonian. Typically you would want to take a peek at the [source code](https://github.com/QuantumKitHub/MPSKitModels.jl/blob/main/src/models/hamiltonians.jl#L379) of `MPSKitModels.jl` to see how these models are defined and tweak it as per your needs. Here we see that applying twisted boundary conditions is equivalent to adding a prefactor of $e^{\pm i\phi}$ in front of the hopping amplitudes. ````julia -function bose_hubbard_model_twisted_bc(elt::Type{<:Number}=ComplexF64, symmetry::Type{<:Sector}=Trivial, - lattice::AbstractLattice=InfiniteChain(1); - cutoff::Integer=5, t=1.0, U=1.0, mu=0.0, phi=0) +function bose_hubbard_model_twisted_bc( + elt::Type{<:Number} = ComplexF64, symmetry::Type{<:Sector} = Trivial, + lattice::AbstractLattice = InfiniteChain(1); + cutoff::Integer = 5, t = 1.0, U = 1.0, mu = 0.0, phi = 0 + ) - a_pm = a_plusmin(elt, symmetry; cutoff=cutoff) - a_mp = a_minplus(elt, symmetry; cutoff=cutoff) - N = a_number(elt, symmetry; cutoff=cutoff) + a_pm = a_plusmin(elt, symmetry; cutoff = cutoff) + a_mp = a_minplus(elt, symmetry; cutoff = cutoff) + N = a_number(elt, symmetry; cutoff = cutoff) interaction_term = N * (N - id(domain(N))) - H = @mpoham begin + return H = @mpoham begin sum(nearest_neighbours(lattice)) do (i, j) - return -t * (exp(1im * phi) * a_pm{i,j} + exp(1im * -phi) * a_mp{i,j}) + return -t * (exp(1im * phi) * a_pm{i, j} + exp(1im * -phi) * a_mp{i, j}) end + - sum(vertices(lattice)) do i + sum(vertices(lattice)) do i return U / 2 * interaction_term{i} - mu * N{i} end end end -function superfluid_stiffness_profile(t, mu, D, cutoff, ϵ=1e-4, npoints=11) +function superfluid_stiffness_profile(t, mu, D, cutoff, ϵ = 1.0e-4, npoints = 11) phis = range(-ϵ, ϵ, npoints) energies = zeros(length(phis)) Threads.@threads for idx in eachindex(phis) - hamiltonian_twisted = bose_hubbard_model_twisted_bc(cutoff=cutoff, t=t, mu=mu, U=1, phi=phis[idx]) + hamiltonian_twisted = bose_hubbard_model_twisted_bc(cutoff = cutoff, t = t, mu = mu, U = 1, phi = phis[idx]) state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) - state_twisted, _, _ = find_groundstate(state, hamiltonian_twisted, VUMPS(; tol=1e-8, verbosity=0)) + state_twisted, _, _ = find_groundstate(state, hamiltonian_twisted, VUMPS(; tol = 1.0e-8, verbosity = 0)) energies[idx] = real(expectation_value(state_twisted, hamiltonian_twisted)) end - plot(phis, energies, lw=2, xlabel="Phase twist per site" * L"(\phi)", ylabel="Ground state energy", title="t = $t | μ = $mu | D = $D | cutoff = $cutoff") + return plot(phis, energies, lw = 2, xlabel = "Phase twist per site" * L"(\phi)", ylabel = "Ground state energy", title = "t = $t | μ = $mu | D = $D | cutoff = $cutoff") end superfluid_stiffness_profile(0.2, 0.3, 5, 4) # superfluid @@ -311,7 +323,7 @@ superfluid_stiffness_profile(0.01, 0.3, 5, 4) # mott insulator ```` ```@raw html - + ``` Now that we know what phases to expect, we can plot the phase diagram by scanning over a range of parameters. In general, one could do better by performing a bisection algorithm for each chemical potential to determine the value of the hopping parameter at the transition point, however the 1D Bose-Hubbard model may have two transition points at the same chemical potential which makes this a bit cumbersome to implement robustly. Furthermore, we stick to using the quasi-condensate density as an order parameter since extracting the superfluid density accurately requires a more robust scheme to compute second derivatives which takes us away from the focus of this tutorial. @@ -321,21 +333,21 @@ cutoff, D = 4, 10 mus = range(0, 0.75, 40) ts = range(0, 0.3, 40) -a_op = a_min(cutoff=cutoff) +a_op = a_min(cutoff = cutoff) order_parameters = zeros(length(ts), length(mus)) Threads.@threads for (i, j) in collect(Iterators.product(eachindex(mus), eachindex(ts))) - hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff=cutoff, U=1, mu=mus[i], t=ts[j]) + hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff = cutoff, U = 1, mu = mus[i], t = ts[j]) init_state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) - state, _, _ = find_groundstate(init_state, hamiltonian, VUMPS(; tol=1e-8, verbosity=0)) + state, _, _ = find_groundstate(init_state, hamiltonian, VUMPS(; tol = 1.0e-8, verbosity = 0)) order_parameters[i, j] = abs(expectation_value(state, 0 => a_op)) end -heatmap(ts, mus, order_parameters, xlabel=L"t/U", ylabel=L"\mu/U", title=L"\langle \hat{a}_i \rangle") +heatmap(ts, mus, order_parameters, xlabel = L"t/U", ylabel = L"\mu/U", title = L"\langle \hat{a}_i \rangle") ```` ```@raw html - + ``` Although the bond dimension here is quite low, we already see the deformation of the Mott insulator lobes to give way to the well known BKT transition that happens at commensurate density. One can go further and estimate the critical exponents using finite-entanglement scaling procedures on the correlation functions, but these may now be performed with ease using what we have learnt in this tutorial. diff --git a/docs/src/examples/quantum1d/8.bose-hubbard/main.ipynb b/docs/src/examples/quantum1d/8.bose-hubbard/main.ipynb index b58889206..e909ec636 100644 --- a/docs/src/examples/quantum1d/8.bose-hubbard/main.ipynb +++ b/docs/src/examples/quantum1d/8.bose-hubbard/main.ipynb @@ -9,7 +9,7 @@ "using Plots, LaTeXStrings\n", "\n", "theme(:wong)\n", - "default(fontfamily=\"Computer Modern\", label=nothing, dpi=100, framestyle=:box)" + "default(fontfamily = \"Computer Modern\", label = nothing, dpi = 100, framestyle = :box)" ], "metadata": {}, "execution_count": null @@ -88,9 +88,9 @@ "outputs": [], "cell_type": "code", "source": [ - "a_op = a_min(cutoff=cutoff) # creates a bosonic annihilation operator without any symmetries\n", + "a_op = a_min(cutoff = cutoff) # creates a bosonic annihilation operator without any symmetries\n", "display(a_op[])\n", - "display((a_op'*a_op)[])" + "display((a_op' * a_op)[])" ], "metadata": {}, "execution_count": null @@ -106,7 +106,7 @@ "outputs": [], "cell_type": "code", "source": [ - "hamiltonian = bose_hubbard_model(InfiniteChain(1), cutoff=cutoff, U=1, mu=0.5, t=0.2) # It is not strictly required to pass InfiniteChain() and is only included for clarity; one may instead pass FiniteChain(N) as well" + "hamiltonian = bose_hubbard_model(InfiniteChain(1), cutoff = cutoff, U = 1, mu = 0.5, t = 0.2) # It is not strictly required to pass InfiniteChain() and is only included for clarity; one may instead pass FiniteChain(N) as well" ], "metadata": {}, "execution_count": null @@ -122,7 +122,7 @@ "outputs": [], "cell_type": "code", "source": [ - "ground_state, _, _ = find_groundstate(initial_state, hamiltonian, VUMPS(tol=1e-6, verbosity=2, maxiter=200))\n", + "ground_state, _, _ = find_groundstate(initial_state, hamiltonian, VUMPS(tol = 1.0e-6, verbosity = 2, maxiter = 200))\n", "println(\"Energy: \", expectation_value(ground_state, hamiltonian))" ], "metadata": {}, @@ -140,14 +140,14 @@ "cell_type": "code", "source": [ "function get_ground_state(mu, t, cutoff, D; kwargs...)\n", - " hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff=cutoff, U=1, mu=mu, t=t)\n", + " hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff = cutoff, U = 1, mu = mu, t = t)\n", " state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D)\n", " state, _, _ = find_groundstate(state, hamiltonian, VUMPS(; kwargs...))\n", "\n", " return state\n", "end\n", "\n", - "ground_state = get_ground_state(0.5, 0.01, cutoff, D; tol=1e-6, verbosity=2, maxiter=500)" + "ground_state = get_ground_state(0.5, 0.01, cutoff, D; tol = 1.0e-6, verbosity = 2, maxiter = 500)" ], "metadata": {}, "execution_count": null @@ -163,8 +163,8 @@ "outputs": [], "cell_type": "code", "source": [ - "plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:50), lw=2, xlabel=\"Site index\", ylabel=\"Correlation function\", yscale=:log10)\n", - "hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls=:dash, c=:black)" + "plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:50), lw = 2, xlabel = \"Site index\", ylabel = \"Correlation function\", yscale = :log10)\n", + "hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls = :dash, c = :black)" ], "metadata": {}, "execution_count": null @@ -180,10 +180,10 @@ "outputs": [], "cell_type": "code", "source": [ - "ground_state = get_ground_state(0.5, 0.2, cutoff, D; tol=1e-6, verbosity=2, maxiter=500)\n", + "ground_state = get_ground_state(0.5, 0.2, cutoff, D; tol = 1.0e-6, verbosity = 2, maxiter = 500)\n", "\n", - "plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:100), lw=2, xlabel=\"Site index\", ylabel=\"Correlation function\", yscale=:log10, xscale=:log10)\n", - "hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls=:dash, c=:black)" + "plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:100), lw = 2, xlabel = \"Site index\", ylabel = \"Correlation function\", yscale = :log10, xscale = :log10)\n", + "hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls = :dash, c = :black)" ], "metadata": {}, "execution_count": null @@ -205,36 +205,42 @@ "states = Vector{InfiniteMPS}(undef, length(Ds))\n", "\n", "Threads.@threads for idx in eachindex(Ds)\n", - " states[idx] = get_ground_state(mu, t, cutoff, Ds[idx], tol=1e-7, verbosity=1, maxiter=500)\n", + " states[idx] = get_ground_state(mu, t, cutoff, Ds[idx], tol = 1.0e-7, verbosity = 1, maxiter = 500)\n", "end\n", "\n", "npoints = 400\n", "two_point_correlation = zeros(length(Ds), npoints)\n", - "a_op = a_min(cutoff=cutoff)\n", + "a_op = a_min(cutoff = cutoff)\n", "\n", "Threads.@threads for idx in eachindex(Ds)\n", " two_point_correlation[idx, :] .= real.(expectation_value(states[idx], (1, i) => a_op' ⊗ a_op) for i in 1:npoints)\n", "end\n", "\n", - "p = plot(framestyle=:box, ylabel=\"Correlation function \" * L\"\\langle a_i^{\\dagger}a_j \\rangle\",\n", - " xlabel=\"Distance \" * L\"|i-j|\", xscale=:log10, yscale=:log10,\n", - " xticks=([10, 100], [\"10\", \"100\"]),\n", - " yticks=([0.25, 0.5, 1.0], [\"0.25\", \"0.5\", \"1.0\"]))\n", - "\n", - "plot!(p, 2:npoints, two_point_correlation[:, 2:end]',\n", - " lab=\"D = \" .* string.(permutedims(Ds)), lw=2)\n", - "\n", - "scatter!(p, Ds, correlation_length.(states),\n", - " ylabel=\"Correlation length\", xlabel=\"Bond dimension\",\n", - " xscale=:log10, yscale=:log10,\n", - " inset=bbox(0.2, 0.51, 0.25, 0.25),\n", - " subplot=2,\n", - " xticks=(20:10:50, string.(20:10:50)),\n", - " yticks=([50, 100], string.([50, 100])),\n", - " xlabelfontsize=8,\n", - " ylabelfontsize=8,\n", - " ylims=[20, 130],\n", - " xlims=[15, 60])" + "p = plot(\n", + " framestyle = :box, ylabel = \"Correlation function \" * L\"\\langle a_i^{\\dagger}a_j \\rangle\",\n", + " xlabel = \"Distance \" * L\"|i-j|\", xscale = :log10, yscale = :log10,\n", + " xticks = ([10, 100], [\"10\", \"100\"]),\n", + " yticks = ([0.25, 0.5, 1.0], [\"0.25\", \"0.5\", \"1.0\"])\n", + ")\n", + "\n", + "plot!(\n", + " p, 2:npoints, two_point_correlation[:, 2:end]',\n", + " lab = \"D = \" .* string.(permutedims(Ds)), lw = 2\n", + ")\n", + "\n", + "scatter!(\n", + " p, Ds, correlation_length.(states),\n", + " ylabel = \"Correlation length\", xlabel = \"Bond dimension\",\n", + " xscale = :log10, yscale = :log10,\n", + " inset = bbox(0.2, 0.51, 0.25, 0.25),\n", + " subplot = 2,\n", + " xticks = (20:10:50, string.(20:10:50)),\n", + " yticks = ([50, 100], string.([50, 100])),\n", + " xlabelfontsize = 8,\n", + " ylabelfontsize = 8,\n", + " ylims = [20, 130],\n", + " xlims = [15, 60]\n", + ")" ], "metadata": {}, "execution_count": null @@ -292,11 +298,15 @@ "ks = range(-0.05, 0.15, 500)\n", "momentum_distribution = map(\n", " ((corr, qc),) -> sum(\n", - " 2 .* cos.(ks' .* (2:npoints)) .* (corr[2:end] .- qc), dims=1) .+ (corr[1] .- qc),\n", - " zip(eachrow(two_point_correlation),\n", - " quasicondensate_density))\n", + " 2 .* cos.(ks' .* (2:npoints)) .* (corr[2:end] .- qc), dims = 1\n", + " ) .+ (corr[1] .- qc),\n", + " zip(\n", + " eachrow(two_point_correlation),\n", + " quasicondensate_density\n", + " )\n", + ")\n", "momentum_distribution = vcat(momentum_distribution...)'\n", - "plot(ks, momentum_distribution, lab=\"D = \" .* string.(permutedims(Ds)), lw=1.5, xlabel=\"Momentum k\", ylabel=L\"\\langle n_k \\rangle\", ylim=[0, 50])" + "plot(ks, momentum_distribution, lab = \"D = \" .* string.(permutedims(Ds)), lw = 1.5, xlabel = \"Momentum k\", ylabel = L\"\\langle n_k \\rangle\", ylim = [0, 50])" ], "metadata": {}, "execution_count": null @@ -318,38 +328,40 @@ "outputs": [], "cell_type": "code", "source": [ - "function bose_hubbard_model_twisted_bc(elt::Type{<:Number}=ComplexF64, symmetry::Type{<:Sector}=Trivial,\n", - " lattice::AbstractLattice=InfiniteChain(1);\n", - " cutoff::Integer=5, t=1.0, U=1.0, mu=0.0, phi=0)\n", + "function bose_hubbard_model_twisted_bc(\n", + " elt::Type{<:Number} = ComplexF64, symmetry::Type{<:Sector} = Trivial,\n", + " lattice::AbstractLattice = InfiniteChain(1);\n", + " cutoff::Integer = 5, t = 1.0, U = 1.0, mu = 0.0, phi = 0\n", + " )\n", "\n", - " a_pm = a_plusmin(elt, symmetry; cutoff=cutoff)\n", - " a_mp = a_minplus(elt, symmetry; cutoff=cutoff)\n", - " N = a_number(elt, symmetry; cutoff=cutoff)\n", + " a_pm = a_plusmin(elt, symmetry; cutoff = cutoff)\n", + " a_mp = a_minplus(elt, symmetry; cutoff = cutoff)\n", + " N = a_number(elt, symmetry; cutoff = cutoff)\n", "\n", " interaction_term = N * (N - id(domain(N)))\n", "\n", - " H = @mpoham begin\n", + " return H = @mpoham begin\n", " sum(nearest_neighbours(lattice)) do (i, j)\n", - " return -t * (exp(1im * phi) * a_pm{i,j} + exp(1im * -phi) * a_mp{i,j})\n", + " return -t * (exp(1im * phi) * a_pm{i, j} + exp(1im * -phi) * a_mp{i, j})\n", " end +\n", - " sum(vertices(lattice)) do i\n", + " sum(vertices(lattice)) do i\n", " return U / 2 * interaction_term{i} - mu * N{i}\n", " end\n", " end\n", "end\n", "\n", - "function superfluid_stiffness_profile(t, mu, D, cutoff, ϵ=1e-4, npoints=11)\n", + "function superfluid_stiffness_profile(t, mu, D, cutoff, ϵ = 1.0e-4, npoints = 11)\n", " phis = range(-ϵ, ϵ, npoints)\n", " energies = zeros(length(phis))\n", "\n", " Threads.@threads for idx in eachindex(phis)\n", - " hamiltonian_twisted = bose_hubbard_model_twisted_bc(cutoff=cutoff, t=t, mu=mu, U=1, phi=phis[idx])\n", + " hamiltonian_twisted = bose_hubbard_model_twisted_bc(cutoff = cutoff, t = t, mu = mu, U = 1, phi = phis[idx])\n", " state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D)\n", - " state_twisted, _, _ = find_groundstate(state, hamiltonian_twisted, VUMPS(; tol=1e-8, verbosity=0))\n", + " state_twisted, _, _ = find_groundstate(state, hamiltonian_twisted, VUMPS(; tol = 1.0e-8, verbosity = 0))\n", " energies[idx] = real(expectation_value(state_twisted, hamiltonian_twisted))\n", " end\n", "\n", - " plot(phis, energies, lw=2, xlabel=\"Phase twist per site\" * L\"(\\phi)\", ylabel=\"Ground state energy\", title=\"t = $t | μ = $mu | D = $D | cutoff = $cutoff\")\n", + " return plot(phis, energies, lw = 2, xlabel = \"Phase twist per site\" * L\"(\\phi)\", ylabel = \"Ground state energy\", title = \"t = $t | μ = $mu | D = $D | cutoff = $cutoff\")\n", "end\n", "\n", "superfluid_stiffness_profile(0.2, 0.3, 5, 4) # superfluid\n", @@ -374,17 +386,17 @@ "mus = range(0, 0.75, 40)\n", "ts = range(0, 0.3, 40)\n", "\n", - "a_op = a_min(cutoff=cutoff)\n", + "a_op = a_min(cutoff = cutoff)\n", "order_parameters = zeros(length(ts), length(mus))\n", "\n", "Threads.@threads for (i, j) in collect(Iterators.product(eachindex(mus), eachindex(ts)))\n", - " hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff=cutoff, U=1, mu=mus[i], t=ts[j])\n", + " hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff = cutoff, U = 1, mu = mus[i], t = ts[j])\n", " init_state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D)\n", - " state, _, _ = find_groundstate(init_state, hamiltonian, VUMPS(; tol=1e-8, verbosity=0))\n", + " state, _, _ = find_groundstate(init_state, hamiltonian, VUMPS(; tol = 1.0e-8, verbosity = 0))\n", " order_parameters[i, j] = abs(expectation_value(state, 0 => a_op))\n", "end\n", "\n", - "heatmap(ts, mus, order_parameters, xlabel=L\"t/U\", ylabel=L\"\\mu/U\", title=L\"\\langle \\hat{a}_i \\rangle\")" + "heatmap(ts, mus, order_parameters, xlabel = L\"t/U\", ylabel = L\"\\mu/U\", title = L\"\\langle \\hat{a}_i \\rangle\")" ], "metadata": {}, "execution_count": null From ec6dd36e48d74685ad2a198764390770a89b4f5d Mon Sep 17 00:00:00 2001 From: leburgel Date: Tue, 18 Nov 2025 11:46:48 +0100 Subject: [PATCH 4/8] Fix some math typesetting, add extra (external) references, reflow paragraphs --- docs/make.jl | 3 +- examples/quantum1d/8.bose-hubbard/main.jl | 161 ++++++++++++++++++---- src/algorithms/toolbox.jl | 2 +- 3 files changed, 137 insertions(+), 29 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index de1a1ed1d..301cb6230 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -31,7 +31,8 @@ links = InterLinks( "TensorOperations" => "https://quantumkithub.github.io/TensorOperations.jl/stable/", "KrylovKit" => "https://jutho.github.io/KrylovKit.jl/stable/", "BlockTensorKit" => "https://lkdvos.github.io/BlockTensorKit.jl/dev/", - "MatrixAlgebraKit" => "https://quantumkithub.github.io/MatrixAlgebraKit.jl/stable/" + "MatrixAlgebraKit" => "https://quantumkithub.github.io/MatrixAlgebraKit.jl/stable/", + "MPSKitModels" => "https://quantumkithub.github.io/MPSKitModels.jl/dev/" ) # include MPSKit in all doctests diff --git a/examples/quantum1d/8.bose-hubbard/main.jl b/examples/quantum1d/8.bose-hubbard/main.jl index 665ac9f2d..f525e9ef4 100644 --- a/examples/quantum1d/8.bose-hubbard/main.jl +++ b/examples/quantum1d/8.bose-hubbard/main.jl @@ -8,11 +8,17 @@ default(fontfamily = "Computer Modern", label = nothing, dpi = 100, framestyle = md""" # 1D Bose-Hubbard model -In this tutorial, we will explore the physics of the one-dimensional Bose–Hubbard model using matrix product states. For the most part, we replicate the results presented in [**Phys. Rev. B 105, 134502**](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.105.134502), which can be consulted for any statements in this tutorial that are not otherwise cited. The Hamiltonian under study is now defined as follows: +In this tutorial, we will explore the physics of the one-dimensional Bose–Hubbard model +using matrix product states. For the most part, we replicate the results presented in +[**Phys. Rev. B 105, +134502**](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.105.134502), which can be +consulted for any statements in this tutorial that are not otherwise cited. The Hamiltonian +under study is now defined as follows: $$H = -t \sum_{i} (\hat{a}_i^{\dagger} \hat{a}_{i+1} + \hat{a}_{i+1}^{\dagger} \hat{a}_i) + \frac{U}{2} \sum_i \hat{n}_i(\hat{n}_i - 1) - \mu \sum_i \hat{n}_i$$ -where the bosonic creation and annihilation operators satisfy the canonical commutation relations (CCR): +where the bosonic creation and annihilation operators satisfy the canonical commutation +relations (CCR): $$[\hat{a}_i, \hat{a}_j^{\dagger}] = \delta_{ij}.$$ @@ -20,7 +26,7 @@ Each lattice site hosts a local Hilbert space corresponding to bosonic occupatio Within this truncated space, the local creation and annihilation operators are represented by finite-dimensional matrices. For example, with cutoff $n_{\text{max}}$, the annihilation operator takes the form -$$ +```math \hat{a} = \begin{bmatrix} 0 & \sqrt{1} & 0 & 0 & \cdots & 0 \\ @@ -30,11 +36,11 @@ $$ 0 & 0 & 0 & \cdots & 0 & \sqrt{n_{\text{max}}} \\ 0 & 0 & 0 & \cdots & 0 & 0 \end{bmatrix} -$$ +``` and the creation operator is simply its Hermitian conjugate, -$$ +```math \hat{a}^\dagger = \begin{bmatrix} 0 & 0 & 0 & \cdots & 0 & 0 \\ @@ -44,22 +50,45 @@ $$ 0 & 0 & 0 & \cdots & 0 & 0 \\ 0 & 0 & 0 & \cdots & \sqrt{n_{\text{max}}} & 0 \end{bmatrix} -$$ +``` The number operator is then given by $$\hat{n} = \hat{a}^\dagger \hat{a} = \mathrm{diag}(0, 1, 2, \ldots, n_{\text{max}}).$$ -Before moving on, notice that the Hamiltonian is uniform and translationally invariant. Typically, such models are studied on a finite chain of $N$ sites with periodic boundary conditions, but this introduces finite-size effects that are rather annoying to deal with. In contrast, the MPS framework allows us to work directly in the thermodynamic limit, avoiding such artifacts. We will follow this line of exploration in this tutorial and leave finite systems for another time. - -In order to work in the thermodynamic limit, we will have to create an `InfiniteMPS`. A complete specification of the MPS requires us to define the physical space and the virtual space of the constituent tensors. At this point is it useful to note that `MPSKit.jl` is powered by `TensorKit.jl` under the hood and has some very generic interfaces in order to allow imposing symmetries of all kinds. As a result, it is sometimes necessary to be a bit more explicit about what we want to do in terms of the vector spaces involved. In this case, we will not consider any symmetries and simply take the most naive approach of working within the `Trivial` sector. The physical space is then `ℂ^(nmax+1)`, and the virtual space is `ℂ^D` where $D$ is some integer chosen to be the bond dimension of the MPS, and `ℂ` is an alias for `ComplexSpace` (typeset as `\bbC`). As $D$ is increased, one increases the amount of entanglement, i.e, quantum correlations that can be captured by the state. +Before moving on, notice that the Hamiltonian is uniform and translationally invariant. +Typically, such models are studied on a finite chain of $N$ sites with periodic boundary +conditions, but this introduces finite-size effects that are rather annoying to deal with. +In contrast, the MPS framework allows us to work directly in the thermodynamic limit, +avoiding such artifacts. We will follow this line of exploration in this tutorial and leave +finite systems for another time. + +In order to work in the thermodynamic limit, we will have to create an +[`InfiniteMPS`](@ref). A complete specification of the MPS requires us to define the +physical space and the virtual space of the constituent tensors. At this point is it useful +to note that `MPSKit.jl` is powered by +[`TensorKit.jl`](https://github.com/QuantumKitHub/TensorKit.jl) under the hood and has some +very generic interfaces in order to allow imposing symmetries of all kinds. As a result, it +is sometimes necessary to be a bit more explicit about what we want to do in terms of the +vector spaces involved. In this case, we will not consider any symmetries and simply take +the most naive approach of working within the [`Trivial`](@extref TensorKitSectors.Trivial) +sector. The physical space is then `ℂ^(nmax+1)`, and the virtual space is `ℂ^D` where $D$ is +some integer chosen to be the bond dimension of the MPS, and `ℂ` is an alias for +[`ComplexSpace`](@extref TensorKit.ComplexSpace) (typeset as `\bbC`). As $D$ is increased, +one increases the amount of entanglement, i.e, quantum correlations that can be captured by +the state. """ cutoff, D = 4, 5 initial_state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) md""" -This simply initializes a tensor filled with random entries (check out the documentation for other useful constructors). Next, we need the creation and annihilation operators. While we could construct them from scratch, here we will use `MPSKitModels.jl` instead that has predefined operators and models for most well-known lattice models. +This simply initializes a tensor filled with random entries (check out the documentation for +other useful constructors). Next, we need the creation and annihilation operators. While we +could construct them from scratch, here we will use +[`MPSKitModels.jl`](https://github.com/QuantumKitHub/MPSKitModels.jl) instead that has +predefined operators and models for most well-known lattice models. In particular, we can +use [`MPSKitModels.a_min`](@extref) to create the bosonic annihilation operator. """ a_op = a_min(cutoff = cutoff) # creates a bosonic annihilation operator without any symmetries @@ -67,20 +96,31 @@ display(a_op[]) display((a_op' * a_op)[]) md""" -The [] accessor lets us see the underlying array and indeed the operators are exactly what we require. Similarly, the Bose Hubbard model is also predefined (although we will construct our own variant later on). + +The [] accessor lets us see the underlying array, and indeed the operators are exactly what +we require. Similarly, the Bose Hubbard model is also predefined in +[`MPSKitModels.bose_hubbard_model`](@extref) also predefined (although we will construct our +own variant later on). + """ hamiltonian = bose_hubbard_model(InfiniteChain(1), cutoff = cutoff, U = 1, mu = 0.5, t = 0.2) # It is not strictly required to pass InfiniteChain() and is only included for clarity; one may instead pass FiniteChain(N) as well md""" -This has created the Hamiltonian operator as a matrix product operator (MPO) which is a convenient form to use in conjunction with MPS. Finally, the ground state optimization may be performed with either `iDMRG` or `VUMPS`. Both should take similar arguments but it is known that VUMPS is typically more efficient for these systems so we proceed with that. +This has created the Hamiltonian operator as a [matrix product operator](@ref +InfiniteMPOHamiltonian) (MPO) which is a convenient form to use in conjunction with MPS. +Finally, the ground state optimization may be performed with either [`iDMRG`](@ref IDMRG) or +[`VUMPS`](@ref). Both should take similar arguments but it is known that VUMPS is typically +more efficient for these systems so we proceed with that. """ ground_state, _, _ = find_groundstate(initial_state, hamiltonian, VUMPS(tol = 1.0e-6, verbosity = 2, maxiter = 200)) println("Energy: ", expectation_value(ground_state, hamiltonian)) md""" -This automatically runs the algorithm until a certain error [measure](https://github.com/QuantumKitHub/MPSKit.jl/blob/main/src/algorithms/toolbox.jl#L42) falls below the specified tolerance or the maximum iterations is reached. Let us wrap all this into a convenient function. +This automatically runs the algorithm until a certain error [measure](@ref +MPSKit.calc_galerkin) falls below the specified tolerance or the maximum iterations is +reached. Let us wrap all this into a convenient function. """ function get_ground_state(mu, t, cutoff, D; kwargs...) @@ -94,14 +134,20 @@ end ground_state = get_ground_state(0.5, 0.01, cutoff, D; tol = 1.0e-6, verbosity = 2, maxiter = 500) md""" -Now that we have the state, we may compute observables using the `expectation_value` function. It typically expects a `Pair`, `(i1, i2, .., ik) => op` where `op` is a `TensorMap` or `InfiniteMPO` acting over `k` sites. In case of the Hamiltonian, it is not necessary to specify the indices as it spans the whole lattice. We can now plot the correlation function $\langle \hat{a}^{\dagger}_i \hat{a}_j\rangle$. + +Now that we have the state, we may compute observables using the [`expectation_value`](@ref) +function. It typically expects a `Pair`, `(i1, i2, .., ik) => op` where `op` is a +`TensorMap` or `InfiniteMPO` acting over `k` sites. In case of the Hamiltonian, it is not +necessary to specify the indices as it spans the whole lattice. We can now plot the +correlation function $\langle \hat{a}^{\dagger}_i \hat{a}_j\rangle$. """ plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:50), lw = 2, xlabel = "Site index", ylabel = "Correlation function", yscale = :log10) hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls = :dash, c = :black) md""" -We see that the correlations drop off exponentially, indicating the existence of a gapped Mott insulating phase. Let us now shift our parameters to probe other phases. +We see that the correlations drop off exponentially, indicating the existence of a gapped +Mott insulating phase. Let us now shift our parameters to probe other phases. """ ground_state = get_ground_state(0.5, 0.2, cutoff, D; tol = 1.0e-6, verbosity = 2, maxiter = 500) @@ -110,7 +156,19 @@ plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls = :dash, c = :black) md""" -In this case, the correlation function drops off algebraically and eventually saturates as $\lim_{i \to \infty}\langle\hat{a}_i^{\dagger} \hat{a}_j\rangle ≈ \langle \hat{a}_i^{\dagger}\rangle \langle \hat{a}_j \rangle = |\langle a_i \rangle|^2 \neq 0$. This is a signature of long-range order and suggests the existence of a Bose-Einstein condensate. However, this is a bit odd since at zero temperature, the Bose Hubbard model is not expected to break any continuous symmetries ($U(1)$ in this case, corresponding to particle number conservation) due to the Mermin-Wagner theorem. The source of this contradiction lies in the fact that the true 1D superfluid ground state is an extended critical phase exhibiting algebraic decay, however, a finite bond-dimension MPS can only capture exponentially decaying correlations. As a result, the finite bond dimension effectively introduces a length scale into the system in a similar manner as finite-size effects. We can see this clearly by increasing the bond dimension. We also see that the correlation length seems to depend algebraically on the bond dimension as expected from finite-entanglement scaling arguments. +In this case, the correlation function drops off algebraically and eventually saturates as +$\lim_{i \to \infty}\langle\hat{a}_i^{\dagger} \hat{a}_j\rangle ≈ \langle \hat{a}_i^{\dagger}\rangle \langle \hat{a}_j \rangle = |\langle a_i \rangle|^2 \neq 0$. +This is a signature of long-range order and suggests the existence of a Bose-Einstein +condensate. However, this is a bit odd since at zero temperature, the Bose Hubbard model is +not expected to break any continuous symmetries ($U(1)$ in this case, corresponding to +particle number conservation) due to the Mermin-Wagner theorem. The source of this +contradiction lies in the fact that the true 1D superfluid ground state is an extended +critical phase exhibiting algebraic decay, however, a finite bond-dimension MPS can only +capture exponentially decaying correlations. As a result, the finite bond dimension +effectively introduces a length scale into the system in a similar manner as finite-size +effects. We can see this clearly by increasing the bond dimension. We also see that the +correlation length seems to depend algebraically on the bond dimension as expected from +finite-entanglement scaling arguments. """ cutoff = 4 @@ -157,7 +215,11 @@ scatter!( ) md""" -This shows that any finite bond dimension MPS necessarily breaks the symmetry of the system, forming a Bose-Einstein condensate which introduces erraneous long-distance behaviour of correlation functions. In case of finite bond dimension, it is thus reasonable to associate the finite expectation value of the field operator to the 'quasicondensate' density of the system which vanishes as $D \to \infty$. +This shows that any finite bond dimension MPS necessarily breaks the symmetry of the system, +forming a Bose-Einstein condensate which introduces erraneous long-distance behaviour of +correlation functions. In case of finite bond dimension, it is thus reasonable to associate +the finite expectation value of the field operator to the 'quasicondensate' density of the +system which vanishes as $D \to \infty$. """ quasicondensate_density = map(state -> abs2(expectation_value(state, (0,) => a_op)), states) @@ -165,14 +227,17 @@ quasicondensate_density = map(state -> abs2(expectation_value(state, (0,) => a_o md""" We may now also visualize the momentum distribution function, which is obtained as the Fourier transform of the single-particle density matrix. Starting from the definition of the momentum occupation operator: -$$\hat{a}_k = \frac{1}{\sqrt{L}} \sum_j e^{-ikj} \hat{a}_j, \qquad +```math +\hat{a}_k = \frac{1}{\sqrt{L}} \sum_j e^{-ikj} \hat{a}_j, \qquad \hat{a}_k^\dagger = \frac{1}{\sqrt{L}} \sum_{j'} e^{ikj'} \hat{a}_{j'}^\dagger -$$ +``` the momentum distribution is -$$\langle \hat{n}_k \rangle = \langle \hat{a}_k^\dagger \hat{a}_k \rangle -= \frac{1}{L} \sum_{j',j} e^{ik(j'-j)} \langle \hat{a}_{j'}^\dagger \hat{a}_j \rangle.$$ +```math +\langle \hat{n}_k \rangle = \langle \hat{a}_k^\dagger \hat{a}_k \rangle += \frac{1}{L} \sum_{j',j} e^{ik(j'-j)} \langle \hat{a}_{j'}^\dagger \hat{a}_j \rangle. +``` For a translationally invariant system, the correlation depends only on the distance $r = j' - j$: @@ -186,7 +251,11 @@ The sum over $j$ yields a factor of $L$, which cancels the prefactor, leading to $$\langle \hat{n}_k \rangle = \sum_{r \in \mathbb{Z}} e^{ikr} \langle \hat{a}_r^\dagger \hat{a}_0 \rangle$$ -However, we know that a finite bond dimension MPS introduces a non-zero quasi-condensate density which would give rise to an $\mathcal{O}(N)$ divergence in the momentum distribution that is not indicative of the true physics of the system. Since we know this contribution vanishes in the infinite bond dimension limit, we instead work with $\langle \hat{a}_r^{\dagger} \hat{a}_0 \rangle_c = \langle \hat{a}_r^{\dagger} \hat{a}_0 \rangle - |\langle \hat{a}\rangle|^2$. +However, we know that a finite bond dimension MPS introduces a non-zero quasi-condensate +density which would give rise to an $\mathcal{O}(N)$ divergence in the momentum distribution +that is not indicative of the true physics of the system. Since we know this contribution +vanishes in the infinite bond dimension limit, we instead work with +$\langle \hat{a}_r^{\dagger} \hat{a}_0 \rangle_c = \langle \hat{a}_r^{\dagger} \hat{a}_0 \rangle - |\langle \hat{a}\rangle|^2$. """ ks = range(-0.05, 0.15, 500) @@ -203,13 +272,40 @@ momentum_distribution = vcat(momentum_distribution...)' plot(ks, momentum_distribution, lab = "D = " .* string.(permutedims(Ds)), lw = 1.5, xlabel = "Momentum k", ylabel = L"\langle n_k \rangle", ylim = [0, 50]) md""" -We see that the density seems to peak around $k=0$, this time seemingly becoming more prominent as $D \to \infty$ which seems to suggest again that there is a condensate. However, going by the Penrose-Onsager criterion, the existence of a condensate can be quantified by requiring the leading eigenvalue of the single particle density matrix (i.e, $\langle \hat{n}_{k=0}\rangle = \sum_j \langle \hat{a}_j^{\dagger} \hat{a}_0\rangle$) to diverge as $O(N)$ in the thermodynamic limit. In this case, since the correlations decay as a power law, there is naturally a divergence at low momenta. But this does not imply the existence of a condensate since the order of divergence is much weaker. However, this does indicate the remnants of some kind of condensation in the 1D model despite the quantum fluctuations, leading to the practical utility of defining the concept of a quasicondensate where there is still a notion of phase coherence over short distances. -What this means for us is that, as far as MPS simulations go, we may still utilize the quasicondensate density as an effective order parameter, although it will be less robust as the bond dimension is increased. Alternatively, we realize that the true phase is characterized as being a superfluid (a concept distinct from Bose-Einstein condensation) and can be identified by a non-zero value of the superfluid stiffness (also known as helicity modulus, $\Upsilon$) as defined by Leggett. Upon applying a phase twist $\Phi$ to the boundaries of the system, a superfluid phase would suffer an increase in energy whereas an insulating phase would not. In the thermodynamic limit, one could show that the boundary conditions may be considered as periodic and instead uniformly distribute the phase across the chain as $\hat{a}_i \to \hat{a}_i e^{i\Phi/L}$. Concretely, in the limit of $\Phi/L \to 0$, we have: +We see that the density seems to peak around $k=0$, this time seemingly becoming more +prominent as $D \to \infty$ which seems to suggest again that there is a condensate. +However, going by the Penrose-Onsager criterion, the existence of a condensate can be +quantified by requiring the leading eigenvalue of the single particle density matrix (i.e, +$\langle \hat{n}_{k=0}\rangle = \sum_j \langle \hat{a}_j^{\dagger} \hat{a}_0\rangle$) to +diverge as $O(N)$ in the thermodynamic limit. In this case, since the correlations decay as +a power law, there is naturally a divergence at low momenta. But this does not imply the +existence of a condensate since the order of divergence is much weaker. However, this does +indicate the remnants of some kind of condensation in the 1D model despite the quantum +fluctuations, leading to the practical utility of defining the concept of a quasicondensate +where there is still a notion of phase coherence over short distances. + +What this means for us is that, as far as MPS simulations go, we may still utilize the +quasicondensate density as an effective order parameter, although it will be less robust as +the bond dimension is increased. Alternatively, we realize that the true phase is +characterized as being a superfluid (a concept distinct from Bose-Einstein condensation) and +can be identified by a non-zero value of the superfluid stiffness (also known as helicity +modulus, $\Upsilon$) as defined by Leggett. Upon applying a phase twist $\Phi$ to the +boundaries of the system, a superfluid phase would suffer an increase in energy whereas an +insulating phase would not. In the thermodynamic limit, one could show that the boundary +conditions may be considered as periodic and instead uniformly distribute the phase across +the chain as $\hat{a}_i \to \hat{a}_i e^{i\Phi/L}$. Concretely, in the limit of +$\Phi/L \to 0$, we have: $$\frac{E[\Phi] - E[0]}{L} \approx \frac{1}{2} \Upsilon(L) \bigg (\frac{\Phi}{L}\bigg)^2 + \cdots$$ -In order to find the ground state under these twisted boundary conditions, we must construct our own variant of the Bose-Hubbard Hamiltonian. Typically you would want to take a peek at the [source code](https://github.com/QuantumKitHub/MPSKitModels.jl/blob/main/src/models/hamiltonians.jl#L379) of `MPSKitModels.jl` to see how these models are defined and tweak it as per your needs. Here we see that applying twisted boundary conditions is equivalent to adding a prefactor of $e^{\pm i\phi}$ in front of the hopping amplitudes. +In order to find the ground state under these twisted boundary conditions, we must construct +our own variant of the Bose-Hubbard Hamiltonian. Typically you would want to take a peek at +the +[source code](https://github.com/QuantumKitHub/MPSKitModels.jl/blob/f4c36d9660a9eab05fa253ffd5c20dc6b7df44cc/src/models/hamiltonians.jl#L379-L409) +of `MPSKitModels.jl` to see how these models are defined and tweak it as per your needs. +Here we see that applying twisted boundary conditions is equivalent to adding a prefactor of +$e^{\pm i\phi}$ in front of the hopping amplitudes. """ function bose_hubbard_model_twisted_bc( @@ -253,7 +349,14 @@ superfluid_stiffness_profile(0.2, 0.3, 5, 4) # superfluid superfluid_stiffness_profile(0.01, 0.3, 5, 4) # mott insulator md""" -Now that we know what phases to expect, we can plot the phase diagram by scanning over a range of parameters. In general, one could do better by performing a bisection algorithm for each chemical potential to determine the value of the hopping parameter at the transition point, however the 1D Bose-Hubbard model may have two transition points at the same chemical potential which makes this a bit cumbersome to implement robustly. Furthermore, we stick to using the quasi-condensate density as an order parameter since extracting the superfluid density accurately requires a more robust scheme to compute second derivatives which takes us away from the focus of this tutorial. +Now that we know what phases to expect, we can plot the phase diagram by scanning over a +range of parameters. In general, one could do better by performing a bisection algorithm for +each chemical potential to determine the value of the hopping parameter at the transition +point, however the 1D Bose-Hubbard model may have two transition points at the same chemical +potential which makes this a bit cumbersome to implement robustly. Furthermore, we stick to +using the quasi-condensate density as an order parameter since extracting the superfluid +density accurately requires a more robust scheme to compute second derivatives which takes +us away from the focus of this tutorial. """ cutoff, D = 4, 10 @@ -273,5 +376,9 @@ end heatmap(ts, mus, order_parameters, xlabel = L"t/U", ylabel = L"\mu/U", title = L"\langle \hat{a}_i \rangle") md""" -Although the bond dimension here is quite low, we already see the deformation of the Mott insulator lobes to give way to the well known BKT transition that happens at commensurate density. One can go further and estimate the critical exponents using finite-entanglement scaling procedures on the correlation functions, but these may now be performed with ease using what we have learnt in this tutorial. +Although the bond dimension here is quite low, we already see the deformation of the Mott +insulator lobes to give way to the well known BKT transition that happens at commensurate +density. One can go further and estimate the critical exponents using finite-entanglement +scaling procedures on the correlation functions, but these may now be performed with ease +using what we have learnt in this tutorial. """ diff --git a/src/algorithms/toolbox.jl b/src/algorithms/toolbox.jl index 150f276fb..bbedda09f 100644 --- a/src/algorithms/toolbox.jl +++ b/src/algorithms/toolbox.jl @@ -43,7 +43,7 @@ Calculate the Galerkin error, which is the error between the solution of the ori Concretely, this is the overlap of the current state with the single-site derivative, projected onto the nullspace of the current state: ```math -\\epsilon = |VL * (VL' * \\frac{above}{\\partial AC_{pos}})| +\\epsilon = \\left|VL ⋅ \\left(VL^{\\dagger} ⋅ \\frac{\\partial \\text{above}}{\\partial AC_{\\text{pos}}}\\right)\\right| ``` """ function calc_galerkin( From 2a01e61593c4dde4ee6432304ca60e8cd49dc46a Mon Sep 17 00:00:00 2001 From: leburgel Date: Tue, 18 Nov 2025 12:52:17 +0100 Subject: [PATCH 5/8] Regerate markdown and notebook --- .../quantum1d/8.bose-hubbard/index.md | 221 +++++++++++++----- .../quantum1d/8.bose-hubbard/main.ipynb | 153 ++++++++++-- examples/Cache.toml | 2 +- 3 files changed, 293 insertions(+), 83 deletions(-) diff --git a/docs/src/examples/quantum1d/8.bose-hubbard/index.md b/docs/src/examples/quantum1d/8.bose-hubbard/index.md index 3396719f3..a0ed93336 100644 --- a/docs/src/examples/quantum1d/8.bose-hubbard/index.md +++ b/docs/src/examples/quantum1d/8.bose-hubbard/index.md @@ -17,11 +17,17 @@ default(fontfamily = "Computer Modern", label = nothing, dpi = 100, framestyle = # 1D Bose-Hubbard model -In this tutorial, we will explore the physics of the one-dimensional Bose–Hubbard model using matrix product states. For the most part, we replicate the results presented in [**Phys. Rev. B 105, 134502**](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.105.134502), which can be consulted for any statements in this tutorial that are not otherwise cited. The Hamiltonian under study is now defined as follows: +In this tutorial, we will explore the physics of the one-dimensional Bose–Hubbard model +using matrix product states. For the most part, we replicate the results presented in +[**Phys. Rev. B 105, +134502**](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.105.134502), which can be +consulted for any statements in this tutorial that are not otherwise cited. The Hamiltonian +under study is now defined as follows: $$H = -t \sum_{i} (\hat{a}_i^{\dagger} \hat{a}_{i+1} + \hat{a}_{i+1}^{\dagger} \hat{a}_i) + \frac{U}{2} \sum_i \hat{n}_i(\hat{n}_i - 1) - \mu \sum_i \hat{n}_i$$ -where the bosonic creation and annihilation operators satisfy the canonical commutation relations (CCR): +where the bosonic creation and annihilation operators satisfy the canonical commutation +relations (CCR): $$[\hat{a}_i, \hat{a}_j^{\dagger}] = \delta_{ij}.$$ @@ -29,7 +35,7 @@ Each lattice site hosts a local Hilbert space corresponding to bosonic occupatio Within this truncated space, the local creation and annihilation operators are represented by finite-dimensional matrices. For example, with cutoff $n_{\text{max}}$, the annihilation operator takes the form -$$ +```math \hat{a} = \begin{bmatrix} 0 & \sqrt{1} & 0 & 0 & \cdots & 0 \\ @@ -39,11 +45,11 @@ $$ 0 & 0 & 0 & \cdots & 0 & \sqrt{n_{\text{max}}} \\ 0 & 0 & 0 & \cdots & 0 & 0 \end{bmatrix} -$$ +``` and the creation operator is simply its Hermitian conjugate, -$$ +```math \hat{a}^\dagger = \begin{bmatrix} 0 & 0 & 0 & \cdots & 0 & 0 \\ @@ -53,15 +59,33 @@ $$ 0 & 0 & 0 & \cdots & 0 & 0 \\ 0 & 0 & 0 & \cdots & \sqrt{n_{\text{max}}} & 0 \end{bmatrix} -$$ +``` The number operator is then given by $$\hat{n} = \hat{a}^\dagger \hat{a} = \mathrm{diag}(0, 1, 2, \ldots, n_{\text{max}}).$$ -Before moving on, notice that the Hamiltonian is uniform and translationally invariant. Typically, such models are studied on a finite chain of $N$ sites with periodic boundary conditions, but this introduces finite-size effects that are rather annoying to deal with. In contrast, the MPS framework allows us to work directly in the thermodynamic limit, avoiding such artifacts. We will follow this line of exploration in this tutorial and leave finite systems for another time. - -In order to work in the thermodynamic limit, we will have to create an `InfiniteMPS`. A complete specification of the MPS requires us to define the physical space and the virtual space of the constituent tensors. At this point is it useful to note that `MPSKit.jl` is powered by `TensorKit.jl` under the hood and has some very generic interfaces in order to allow imposing symmetries of all kinds. As a result, it is sometimes necessary to be a bit more explicit about what we want to do in terms of the vector spaces involved. In this case, we will not consider any symmetries and simply take the most naive approach of working within the `Trivial` sector. The physical space is then `ℂ^(nmax+1)`, and the virtual space is `ℂ^D` where $D$ is some integer chosen to be the bond dimension of the MPS, and `ℂ` is an alias for `ComplexSpace` (typeset as `\bbC`). As $D$ is increased, one increases the amount of entanglement, i.e, quantum correlations that can be captured by the state. +Before moving on, notice that the Hamiltonian is uniform and translationally invariant. +Typically, such models are studied on a finite chain of $N$ sites with periodic boundary +conditions, but this introduces finite-size effects that are rather annoying to deal with. +In contrast, the MPS framework allows us to work directly in the thermodynamic limit, +avoiding such artifacts. We will follow this line of exploration in this tutorial and leave +finite systems for another time. + +In order to work in the thermodynamic limit, we will have to create an +[`InfiniteMPS`](@ref). A complete specification of the MPS requires us to define the +physical space and the virtual space of the constituent tensors. At this point is it useful +to note that `MPSKit.jl` is powered by +[`TensorKit.jl`](https://github.com/QuantumKitHub/TensorKit.jl) under the hood and has some +very generic interfaces in order to allow imposing symmetries of all kinds. As a result, it +is sometimes necessary to be a bit more explicit about what we want to do in terms of the +vector spaces involved. In this case, we will not consider any symmetries and simply take +the most naive approach of working within the [`Trivial`](@extref TensorKitSectors.Trivial) +sector. The physical space is then `ℂ^(nmax+1)`, and the virtual space is `ℂ^D` where $D$ is +some integer chosen to be the bond dimension of the MPS, and `ℂ` is an alias for +[`ComplexSpace`](@extref TensorKit.ComplexSpace) (typeset as `\bbC`). As $D$ is increased, +one increases the amount of entanglement, i.e, quantum correlations that can be captured by +the state. ````julia cutoff, D = 4, 5 @@ -69,15 +93,21 @@ initial_state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) ```` ```` -single site InfiniteMPS: -│ ⋮ -│ C[1]: TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 1, Vector{ComplexF64}}(ComplexF64[0.980773+0.0im, 0.146227-0.0278782im, 0.105333-0.010143im, 0.0620678+0.00292624im, 0.0241261-0.00222812im, 0.0+0.0im, 0.0101049+0.0im, 0.00743715-0.00597138im, -0.0011665+0.00106208im, -0.00256436-0.000271479im, 0.0+0.0im, 0.0+0.0im, 0.00580994+0.0im, -0.00381868-0.00101898im, -0.00128419+0.00100647im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.00397274+0.0im, 0.00234497-0.000379053im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.00198522+0.0im], ℂ^5 ← ℂ^5) -├── AL[1]: TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[0.375885+0.305768im, 0.0768109+0.0330615im, 0.0460733+0.0360179im, 0.0240341+0.025383im, 0.00804042+0.0057914im, 0.213348+0.271799im, 0.0517906+0.0297947im, 0.0138421+0.0188126im, 0.00857498+0.0239495im, 0.00590194+0.0070995im, 0.393411+0.386699im, 0.0781762+0.0467087im, 0.0547209+0.0259146im, 0.0239884+0.026463im, 0.0118271+0.00963569im, 0.242906+0.273226im, 0.0539362+0.037471im, 0.0475866+0.0266459im, 0.0156312+0.0217225im, 0.0105492+0.00484183im, 0.195707+0.347509im, 0.0485176+0.057361im, 0.0399442+0.0476088im, 0.0108898+0.0292134im, 0.00742594+0.00871371im, -0.405789-0.094101im, -0.00619161+0.0455026im, -0.0225547-0.107661im, -0.0576217-0.0382695im, -0.0256341-0.0202461im, 0.480858+0.520426im, 0.000291008+0.091172im, 0.139655+0.119484im, 0.0651289-0.0131326im, 0.017084+0.0264884im, -0.21137-0.149131im, -0.0469705-0.104014im, -0.154879-0.00320825im, -0.0192648-0.007215im, -0.0348013-0.0153511im, 0.0504188+0.0087985im, -0.0921246-0.012542im, -0.120437+0.0817629im, 0.0133867-0.0348062im, -0.0271514-0.0101404im, 0.294447+0.113036im, 0.0607909-0.0783091im, -0.0267377-0.0718419im, 0.0296402-0.0108425im, -0.00573401-0.0122223im, 0.0103164-0.142274im, -0.0646881-0.0987699im, -0.0188412+0.0287894im, 0.0265472-0.0232285im, -0.00197098+0.00516289im, 0.0692401-0.176396im, 0.0426622-0.0273936im, -0.00973839-0.0308723im, 0.0120297+0.0305825im, 0.0208432-0.00767607im, -0.251962-0.320331im, -0.0961427-0.00484129im, 0.0287567-0.0230259im, 0.0149031-0.0233169im, 0.0280275-0.00709338im, 0.371239+0.719013im, 0.114536+0.131076im, 0.123631+0.02352im, -0.0317699+0.0460085im, 0.00597557+0.0348219im, -0.125585+0.0936419im, -0.0818155+0.0358405im, -0.0724073+0.022764im, -0.0116091-0.00616551im, 0.00110864+0.00807954im, 0.40961+0.365138im, -0.128697+0.130384im, 0.118639+0.107234im, 0.0506963-0.0355463im, 0.0413307+0.0240077im, -0.125027+0.164525im, -0.113642-0.0873257im, -0.137262-0.00708175im, 0.0456722+0.0461297im, 0.00611336-0.0213785im, -0.39764-0.209872im, -0.00809221+0.0215245im, -0.026733-0.00806058im, -0.0233946-0.0571397im, 0.0202147-0.0273948im, -0.0566267+0.0987188im, -0.0429491-0.0913588im, -0.179772-0.115045im, 0.0264671+0.0644071im, 0.0164148+0.00709543im, 0.353547-0.317267im, 0.108454-0.00632667im, 0.114669-0.10561im, 0.0304626+0.0034387im, 0.041428-0.00870923im, -0.113545-0.159728im, 0.0926693-0.0340042im, 0.006828-0.0369635im, -0.121696+0.0542231im, -0.0053801+0.0419627im, -0.398247-0.168103im, -0.14182+0.206577im, 0.0674782+0.207936im, -0.0729256-0.175661im, -0.0606546+0.0189888im, 0.0538273+0.214954im, -0.183298-0.00539617im, -0.1011+0.0842453im, -0.052954+0.0404571im, -0.0774749+0.0641526im, -0.0490675+0.214829im, -0.00844192+0.151745im, 0.105707+0.130237im, 0.0499562-0.0486145im, 0.0121478-0.0154518im, 0.572407-0.0611613im, -0.0212294-0.128364im, -0.00861733-0.109909im, -0.00619022-0.07167im, -0.0857339-0.0245475im], (ℂ^5 ⊗ ℂ^5) ← ℂ^5) -│ ⋮ +1-site InfiniteMPS(ComplexF64, TensorKit.ComplexSpace) with maximal dimension 5: +| ⋮ +| ℂ^5 +├─[1]─ ℂ^5 +│ ℂ^5 +| ⋮ ```` -This simply initializes a tensor filled with random entries (check out the documentation for other useful constructors). Next, we need the creation and annihilation operators. While we could construct them from scratch, here we will use `MPSKitModels.jl` instead that has predefined operators and models for most well-known lattice models. +This simply initializes a tensor filled with random entries (check out the documentation for +other useful constructors). Next, we need the creation and annihilation operators. While we +could construct them from scratch, here we will use +[`MPSKitModels.jl`](https://github.com/QuantumKitHub/MPSKitModels.jl) instead that has +predefined operators and models for most well-known lattice models. In particular, we can +use [`MPSKitModels.a_min`](@extref) to create the bosonic annihilation operator. ````julia a_op = a_min(cutoff = cutoff) # creates a bosonic annihilation operator without any symmetries @@ -85,21 +115,30 @@ display(a_op[]) display((a_op' * a_op)[]) ```` -The [] accessor lets us see the underlying array and indeed the operators are exactly what we require. Similarly, the Bose Hubbard model is also predefined (although we will construct our own variant later on). +The [] accessor lets us see the underlying array, and indeed the operators are exactly what +we require. Similarly, the Bose Hubbard model is also predefined in +[`MPSKitModels.bose_hubbard_model`](@extref) also predefined (although we will construct our +own variant later on). ````julia hamiltonian = bose_hubbard_model(InfiniteChain(1), cutoff = cutoff, U = 1, mu = 0.5, t = 0.2) # It is not strictly required to pass InfiniteChain() and is only included for clarity; one may instead pass FiniteChain(N) as well ```` ```` -single site InfiniteMPOHamiltonian{MPSKit.JordanMPOTensor{ComplexF64, TensorKit.ComplexSpace, Union{TensorKit.BraidingTensor{ComplexF64, TensorKit.ComplexSpace}, TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 2, Vector{ComplexF64}}}, TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}, TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 2, Vector{ComplexF64}}, TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 1, Vector{ComplexF64}}}}: -╷ ⋮ -┼ W[1]: JordanMPOTensor{ComplexF64, ComplexSpace, Union{BraidingTensor{ComplexF64, ComplexSpace}, TensorMap{ComplexF64, ComplexSpace, 2, 2, Vector{ComplexF64}}}, TensorMap{ComplexF64, ComplexSpace, 2, 1, Vector{ComplexF64}}, TensorMap{ComplexF64, ComplexSpace, 1, 2, Vector{ComplexF64}}, TensorMap{ComplexF64, ComplexSpace, 1, 1, Vector{ComplexF64}}}(((ℂ^1 ⊞ ℂ^2 ⊞ ℂ^1) ⊗ ⊞(ℂ^5)) ← (⊞(ℂ^5) ⊗ (ℂ^1 ⊞ ℂ^2 ⊞ ℂ^1)), SparseBlockTensorMap{Union{BraidingTensor{ComplexF64, ComplexSpace}, TensorMap{ComplexF64, ComplexSpace, 2, 2, Vector{ComplexF64}}}, ComplexF64, ComplexSpace, 2, 2, 4}(Dict{CartesianIndex{4}, Union{BraidingTensor{ComplexF64, ComplexSpace}, TensorMap{ComplexF64, ComplexSpace, 2, 2, Vector{ComplexF64}}}}(), (⊞(ℂ^2) ⊗ ⊞(ℂ^5)) ← (⊞(ℂ^5) ⊗ ⊞(ℂ^2))), SparseBlockTensorMap{TensorMap{ComplexF64, ComplexSpace, 2, 1, Vector{ComplexF64}}, ComplexF64, ComplexSpace, 2, 1, 3}(Dict{CartesianIndex{3}, TensorMap{ComplexF64, ComplexSpace, 2, 1, Vector{ComplexF64}}}(CartesianIndex(1, 1, 1)=>TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 3.16228+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 3.16228+0.0im, 9.66041e-16+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 4.47214+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, -5.54668e-31+0.0im, 4.47214+0.0im, -2.08077e-16+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 5.47723+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 5.47723+0.0im, -7.40208e-16+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 6.32456+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 6.32456+0.0im, 3.05151e-16+0.0im, 0.0+0.0im, 0.0+0.0im], (ℂ^2 ⊗ ℂ^5) ← ℂ^5)), (⊞(ℂ^2) ⊗ ⊞(ℂ^5)) ← ⊞(ℂ^5)), SparseBlockTensorMap{TensorMap{ComplexF64, ComplexSpace, 1, 2, Vector{ComplexF64}}, ComplexF64, ComplexSpace, 1, 2, 3}(Dict{CartesianIndex{3}, TensorMap{ComplexF64, ComplexSpace, 1, 2, Vector{ComplexF64}}}(CartesianIndex(1, 1, 1)=>TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 2, Vector{ComplexF64}}(ComplexF64[0.0+0.0im, -0.0632456+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, -0.0894427+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, -0.109545+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, -0.126491+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, -3.40958e-18+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, -0.0632456+0.0im, 0.0+0.0im, -3.19295e-18+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, -0.0894427+0.0im, 0.0+0.0im, -1.11022e-17+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, -0.109545+0.0im, 0.0+0.0im, 5.55112e-18+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, -0.126491+0.0im, 0.0+0.0im], ℂ^5 ← (ℂ^5 ⊗ ℂ^2))), ⊞(ℂ^5) ← (⊞(ℂ^5) ⊗ ⊞(ℂ^2))), SparseBlockTensorMap{TensorMap{ComplexF64, ComplexSpace, 1, 1, Vector{ComplexF64}}, ComplexF64, ComplexSpace, 1, 1, 2}(Dict{CartesianIndex{2}, TensorMap{ComplexF64, ComplexSpace, 1, 1, Vector{ComplexF64}}}(CartesianIndex(1, 1)=>TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 1, Vector{ComplexF64}}(ComplexF64[0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, -0.5+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 1.5+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 4.0+0.0im], ℂ^5 ← ℂ^5)), ⊞(ℂ^5) ← ⊞(ℂ^5))) -╵ ⋮ +1-site InfiniteMPOHamiltonian(ComplexF64, TensorKit.ComplexSpace) with maximal dimension 4: +| ⋮ +| (ℂ^1 ⊞ ℂ^2 ⊞ ℂ^1) +┼─[1]─ ℂ^5 +│ (ℂ^1 ⊞ ℂ^2 ⊞ ℂ^1) +| ⋮ ```` -This has created the Hamiltonian operator as a matrix product operator (MPO) which is a convenient form to use in conjunction with MPS. Finally, the ground state optimization may be performed with either `iDMRG` or `VUMPS`. Both should take similar arguments but it is known that VUMPS is typically more efficient for these systems so we proceed with that. +This has created the Hamiltonian operator as a [matrix product operator](@ref +InfiniteMPOHamiltonian) (MPO) which is a convenient form to use in conjunction with MPS. +Finally, the ground state optimization may be performed with either [`iDMRG`](@ref IDMRG) or +[`VUMPS`](@ref). Both should take similar arguments but it is known that VUMPS is typically +more efficient for these systems so we proceed with that. ````julia ground_state, _, _ = find_groundstate(initial_state, hamiltonian, VUMPS(tol = 1.0e-6, verbosity = 2, maxiter = 200)) @@ -107,13 +146,15 @@ println("Energy: ", expectation_value(ground_state, hamiltonian)) ```` ```` -[ Info: VUMPS init: obj = +3.877905708364e-01 err = 6.1967e-01 -[ Info: VUMPS conv 33: obj = -6.757777651162e-01 err = 8.3203601611e-07 time = 9.16 sec -Energy: -0.675777765116157 + 6.114900252818245e-17im +[ Info: VUMPS init: obj = +7.210229722103e-01 err = 5.9367e-01 +[ Info: VUMPS conv 47: obj = -6.757777651163e-01 err = 8.0696299728e-07 time = 0.25 sec +Energy: -0.6757777651163069 + 9.5509045594669e-18im ```` -This automatically runs the algorithm until a certain error [measure](https://github.com/QuantumKitHub/MPSKit.jl/blob/main/src/algorithms/toolbox.jl#L42) falls below the specified tolerance or the maximum iterations is reached. Let us wrap all this into a convenient function. +This automatically runs the algorithm until a certain error [measure](@ref +MPSKit.calc_galerkin) falls below the specified tolerance or the maximum iterations is +reached. Let us wrap all this into a convenient function. ````julia function get_ground_state(mu, t, cutoff, D; kwargs...) @@ -128,15 +169,20 @@ ground_state = get_ground_state(0.5, 0.01, cutoff, D; tol = 1.0e-6, verbosity = ```` ```` -single site InfiniteMPS: -│ ⋮ -│ C[1]: TensorMap{ComplexF64, TensorKit.ComplexSpace, 1, 1, Vector{ComplexF64}}(ComplexF64[0.259222+0.0im, 0.246018-0.345952im, 0.610738+0.473339im, -0.148954-0.282818im, -0.211517+0.0836342im, 0.0+0.0im, 0.0136571+0.0im, -0.00924763+0.0239744im, 0.0102051-0.00219558im, 0.00373234-0.000401849im, 0.0+0.0im, 0.0+0.0im, 0.0227273+0.0im, -0.00604063-0.00171163im, 0.000546388+0.0105763im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 3.50256e-6+0.0im, 2.80657e-6-1.34091e-6im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 0.0+0.0im, 4.26623e-6+0.0im], ℂ^5 ← ℂ^5) -├── AL[1]: TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 1, Vector{ComplexF64}}(ComplexF64[-0.087641+0.413009im, 0.214853+0.0747518im, 0.405403+0.3661im, 0.104705-0.379309im, 0.00466853-0.282206im, 0.0435006+0.0621785im, 0.122043-0.0029551im, 0.0130668+0.192364im, 0.0341526-0.0779833im, -0.0402589-0.0307335im, 0.0439795+0.0220837im, 0.077068-0.0405486im, 0.0546722+0.118445im, 0.00333692-0.0560385im, -0.0295748-0.00984742im, 0.0631303-0.0715436im, -0.0365518-0.15324im, 0.283061-0.0553091im, -0.115614-0.0282307im, -0.0282239+0.0789666im, 0.00187513-0.00134459im, 0.000525481-0.000776215im, -0.00180619-3.04145e-5im, -0.000788208+0.000901196im, -0.00130384+0.00111009im, 0.237015-0.0571072im, 0.166803-0.0622261im, -0.197895+0.172091im, -0.0499605+0.0113682im, -0.192109+0.0651415im, -0.021684+0.115554im, 0.121816+0.150109im, -0.260216+0.190125im, 0.126714-0.0298417im, -0.0158932-0.0920934im, 0.128831+0.120805im, 0.271616-0.0633097im, 0.102193+0.525133im, 0.0494367-0.218492im, -0.157755-0.0575024im, -0.0885049+0.036656im, -0.0422998+0.154893im, -0.276699-0.0966281im, 0.0859221+0.0802997im, 0.061836-0.0591176im, -0.00232231+0.00025207im, -0.000850451+0.000423277im, 0.00157813+0.000926226im, 0.0011447-0.000402135im, 0.0017028-0.000329999im, -0.19475+0.104128im, -0.127458+0.218784im, -0.197587-0.269976im, 0.0936674+0.0793825im, 0.154231-0.116272im, 0.183613+0.0437624im, 0.241292-0.214837im, 0.38749+0.468974im, -0.058015-0.238825im, -0.164866+0.0298767im, -0.012259+0.0545234im, 0.0492084+0.0688762im, -0.140159+0.0922497im, 0.0632274-0.0161014im, -0.014811-0.0562813im, -0.071732-0.0151331im, -0.089349+0.0791879im, -0.134631-0.168392im, 0.0228751+0.0856916im, 0.0637356-0.0103791im, -0.00155239-0.000769454im, -0.000700543-7.51919e-5im, 0.000616444+0.0012093im, 0.000875375+0.00020615im, 0.00119626+0.000473548im, -0.210309+0.314047im, -0.146014+0.176102im, 0.267523-0.412974im, -0.0061229-0.115727im, 0.153847-0.254992im, -0.0890551+0.0141788im, -0.0472752+0.120377im, -0.215729-0.106406im, 0.0743765+0.0813334im, 0.0811442-0.0391826im, 0.0468722-0.00895011im, 0.0236013-0.0654369im, 0.152916+0.0731364im, -0.0483316-0.0501106im, -0.050022+0.038829im, -0.0458048+0.137925im, 0.13997+0.192221im, -0.360059+0.239392im, 0.176409-0.0288231im, -0.00697949-0.127442im, -0.00182017+0.00294987im, -0.000279671+0.00138125im, 0.00252286-0.00100342im, 0.000572153-0.00170456im, 0.00116508-0.00229521im, -0.123356-0.0834066im, -0.205016+0.0435919im, -0.0418691-0.332833im, -0.0181937+0.139342im, 0.12447+0.0348314im, -0.0212312-0.0617839im, -0.0831529-0.0234744im, 0.0636293-0.159741im, -0.0468686+0.0599996im, 0.0460042+0.0572296im, 0.158391+0.0558601im, 0.241495-0.160082im, 0.252877+0.417002im, -0.0178357-0.198992im, -0.12711-0.000907711im, 0.11528-0.0948115im, -0.0221029-0.24102im, 0.439961-0.0282941im, -0.17246-0.066882im, -0.0626087+0.113866im, 0.00313034-0.00165794im, 0.000960262-0.00107186im, -0.00274142-0.000420511im, -0.0013851+0.00120634im, -0.0022133+0.00141929im], (ℂ^5 ⊗ ℂ^5) ← ℂ^5) -│ ⋮ +1-site InfiniteMPS(ComplexF64, TensorKit.ComplexSpace) with maximal dimension 5: +| ⋮ +| ℂ^5 +├─[1]─ ℂ^5 +│ ℂ^5 +| ⋮ ```` -Now that we have the state, we may compute observables using the `expectation_value` function. It typically expects a `Pair`, `(i1, i2, .., ik) => op` where `op` is a `TensorMap` or `InfiniteMPO` acting over `k` sites. In case of the Hamiltonian, it is not necessary to specify the indices as it spans the whole lattice. We can now plot the correlation function $\langle \hat{a}^{\dagger}_i \hat{a}_j\rangle$. +Now that we have the state, we may compute observables using the [`expectation_value`](@ref) +function. It typically expects a `Pair`, `(i1, i2, .., ik) => op` where `op` is a +`TensorMap` or `InfiniteMPO` acting over `k` sites. In case of the Hamiltonian, it is not +necessary to specify the indices as it spans the whole lattice. We can now plot the +correlation function $\langle \hat{a}^{\dagger}_i \hat{a}_j\rangle$. ````julia plot(map(i -> real.(expectation_value(ground_state, (0, i) => a_op' ⊗ a_op)), 1:50), lw = 2, xlabel = "Site index", ylabel = "Correlation function", yscale = :log10) @@ -144,10 +190,11 @@ hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls = :dash, c = :b ```` ```@raw html - + ``` -We see that the correlations drop off exponentially, indicating the existence of a gapped Mott insulating phase. Let us now shift our parameters to probe other phases. +We see that the correlations drop off exponentially, indicating the existence of a gapped +Mott insulating phase. Let us now shift our parameters to probe other phases. ````julia ground_state = get_ground_state(0.5, 0.2, cutoff, D; tol = 1.0e-6, verbosity = 2, maxiter = 500) @@ -160,7 +207,19 @@ hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls = :dash, c = :b ``` -In this case, the correlation function drops off algebraically and eventually saturates as $\lim_{i \to \infty}\langle\hat{a}_i^{\dagger} \hat{a}_j\rangle ≈ \langle \hat{a}_i^{\dagger}\rangle \langle \hat{a}_j \rangle = |\langle a_i \rangle|^2 \neq 0$. This is a signature of long-range order and suggests the existence of a Bose-Einstein condensate. However, this is a bit odd since at zero temperature, the Bose Hubbard model is not expected to break any continuous symmetries ($U(1)$ in this case, corresponding to particle number conservation) due to the Mermin-Wagner theorem. The source of this contradiction lies in the fact that the true 1D superfluid ground state is an extended critical phase exhibiting algebraic decay, however, a finite bond-dimension MPS can only capture exponentially decaying correlations. As a result, the finite bond dimension effectively introduces a length scale into the system in a similar manner as finite-size effects. We can see this clearly by increasing the bond dimension. We also see that the correlation length seems to depend algebraically on the bond dimension as expected from finite-entanglement scaling arguments. +In this case, the correlation function drops off algebraically and eventually saturates as +$\lim_{i \to \infty}\langle\hat{a}_i^{\dagger} \hat{a}_j\rangle ≈ \langle \hat{a}_i^{\dagger}\rangle \langle \hat{a}_j \rangle = |\langle a_i \rangle|^2 \neq 0$. +This is a signature of long-range order and suggests the existence of a Bose-Einstein +condensate. However, this is a bit odd since at zero temperature, the Bose Hubbard model is +not expected to break any continuous symmetries ($U(1)$ in this case, corresponding to +particle number conservation) due to the Mermin-Wagner theorem. The source of this +contradiction lies in the fact that the true 1D superfluid ground state is an extended +critical phase exhibiting algebraic decay, however, a finite bond-dimension MPS can only +capture exponentially decaying correlations. As a result, the finite bond dimension +effectively introduces a length scale into the system in a similar manner as finite-size +effects. We can see this clearly by increasing the bond dimension. We also see that the +correlation length seems to depend algebraically on the bond dimension as expected from +finite-entanglement scaling arguments. ````julia cutoff = 4 @@ -208,10 +267,14 @@ scatter!( ```` ```@raw html - + ``` -This shows that any finite bond dimension MPS necessarily breaks the symmetry of the system, forming a Bose-Einstein condensate which introduces erraneous long-distance behaviour of correlation functions. In case of finite bond dimension, it is thus reasonable to associate the finite expectation value of the field operator to the 'quasicondensate' density of the system which vanishes as $D \to \infty$. +This shows that any finite bond dimension MPS necessarily breaks the symmetry of the system, +forming a Bose-Einstein condensate which introduces erraneous long-distance behaviour of +correlation functions. In case of finite bond dimension, it is thus reasonable to associate +the finite expectation value of the field operator to the 'quasicondensate' density of the +system which vanishes as $D \to \infty$. ````julia quasicondensate_density = map(state -> abs2(expectation_value(state, (0,) => a_op)), states) @@ -219,25 +282,28 @@ quasicondensate_density = map(state -> abs2(expectation_value(state, (0,) => a_o ```` 7-element Vector{Float64}: - 0.30974277582170107 - 0.28814776444605805 - 0.27020016270456915 - 0.257127287513971 - 0.24685386110383362 - 0.23539802368529517 - 0.22799655546142042 + 0.3097427747270374 + 0.28814773216310047 + 0.2701999836479986 + 0.25712727854071404 + 0.24685385462377543 + 0.23539744783500957 + 0.2279965610349862 ```` We may now also visualize the momentum distribution function, which is obtained as the Fourier transform of the single-particle density matrix. Starting from the definition of the momentum occupation operator: -$$\hat{a}_k = \frac{1}{\sqrt{L}} \sum_j e^{-ikj} \hat{a}_j, \qquad +```math +\hat{a}_k = \frac{1}{\sqrt{L}} \sum_j e^{-ikj} \hat{a}_j, \qquad \hat{a}_k^\dagger = \frac{1}{\sqrt{L}} \sum_{j'} e^{ikj'} \hat{a}_{j'}^\dagger -$$ +``` the momentum distribution is -$$\langle \hat{n}_k \rangle = \langle \hat{a}_k^\dagger \hat{a}_k \rangle -= \frac{1}{L} \sum_{j',j} e^{ik(j'-j)} \langle \hat{a}_{j'}^\dagger \hat{a}_j \rangle.$$ +```math +\langle \hat{n}_k \rangle = \langle \hat{a}_k^\dagger \hat{a}_k \rangle += \frac{1}{L} \sum_{j',j} e^{ik(j'-j)} \langle \hat{a}_{j'}^\dagger \hat{a}_j \rangle. +``` For a translationally invariant system, the correlation depends only on the distance $r = j' - j$: @@ -251,7 +317,11 @@ The sum over $j$ yields a factor of $L$, which cancels the prefactor, leading to $$\langle \hat{n}_k \rangle = \sum_{r \in \mathbb{Z}} e^{ikr} \langle \hat{a}_r^\dagger \hat{a}_0 \rangle$$ -However, we know that a finite bond dimension MPS introduces a non-zero quasi-condensate density which would give rise to an $\mathcal{O}(N)$ divergence in the momentum distribution that is not indicative of the true physics of the system. Since we know this contribution vanishes in the infinite bond dimension limit, we instead work with $\langle \hat{a}_r^{\dagger} \hat{a}_0 \rangle_c = \langle \hat{a}_r^{\dagger} \hat{a}_0 \rangle - |\langle \hat{a}\rangle|^2$. +However, we know that a finite bond dimension MPS introduces a non-zero quasi-condensate +density which would give rise to an $\mathcal{O}(N)$ divergence in the momentum distribution +that is not indicative of the true physics of the system. Since we know this contribution +vanishes in the infinite bond dimension limit, we instead work with +$\langle \hat{a}_r^{\dagger} \hat{a}_0 \rangle_c = \langle \hat{a}_r^{\dagger} \hat{a}_0 \rangle - |\langle \hat{a}\rangle|^2$. ````julia ks = range(-0.05, 0.15, 500) @@ -269,16 +339,42 @@ plot(ks, momentum_distribution, lab = "D = " .* string.(permutedims(Ds)), lw = 1 ```` ```@raw html - + ``` -We see that the density seems to peak around $k=0$, this time seemingly becoming more prominent as $D \to \infty$ which seems to suggest again that there is a condensate. However, going by the Penrose-Onsager criterion, the existence of a condensate can be quantified by requiring the leading eigenvalue of the single particle density matrix (i.e, $\langle \hat{n}_{k=0}\rangle = \sum_j \langle \hat{a}_j^{\dagger} \hat{a}_0\rangle$) to diverge as $O(N)$ in the thermodynamic limit. In this case, since the correlations decay as a power law, there is naturally a divergence at low momenta. But this does not imply the existence of a condensate since the order of divergence is much weaker. However, this does indicate the remnants of some kind of condensation in the 1D model despite the quantum fluctuations, leading to the practical utility of defining the concept of a quasicondensate where there is still a notion of phase coherence over short distances. - -What this means for us is that, as far as MPS simulations go, we may still utilize the quasicondensate density as an effective order parameter, although it will be less robust as the bond dimension is increased. Alternatively, we realize that the true phase is characterized as being a superfluid (a concept distinct from Bose-Einstein condensation) and can be identified by a non-zero value of the superfluid stiffness (also known as helicity modulus, $\Upsilon$) as defined by Leggett. Upon applying a phase twist $\Phi$ to the boundaries of the system, a superfluid phase would suffer an increase in energy whereas an insulating phase would not. In the thermodynamic limit, one could show that the boundary conditions may be considered as periodic and instead uniformly distribute the phase across the chain as $\hat{a}_i \to \hat{a}_i e^{i\Phi/L}$. Concretely, in the limit of $\Phi/L \to 0$, we have: +We see that the density seems to peak around $k=0$, this time seemingly becoming more +prominent as $D \to \infty$ which seems to suggest again that there is a condensate. +However, going by the Penrose-Onsager criterion, the existence of a condensate can be +quantified by requiring the leading eigenvalue of the single particle density matrix (i.e, +$\langle \hat{n}_{k=0}\rangle = \sum_j \langle \hat{a}_j^{\dagger} \hat{a}_0\rangle$) to +diverge as $O(N)$ in the thermodynamic limit. In this case, since the correlations decay as +a power law, there is naturally a divergence at low momenta. But this does not imply the +existence of a condensate since the order of divergence is much weaker. However, this does +indicate the remnants of some kind of condensation in the 1D model despite the quantum +fluctuations, leading to the practical utility of defining the concept of a quasicondensate +where there is still a notion of phase coherence over short distances. + +What this means for us is that, as far as MPS simulations go, we may still utilize the +quasicondensate density as an effective order parameter, although it will be less robust as +the bond dimension is increased. Alternatively, we realize that the true phase is +characterized as being a superfluid (a concept distinct from Bose-Einstein condensation) and +can be identified by a non-zero value of the superfluid stiffness (also known as helicity +modulus, $\Upsilon$) as defined by Leggett. Upon applying a phase twist $\Phi$ to the +boundaries of the system, a superfluid phase would suffer an increase in energy whereas an +insulating phase would not. In the thermodynamic limit, one could show that the boundary +conditions may be considered as periodic and instead uniformly distribute the phase across +the chain as $\hat{a}_i \to \hat{a}_i e^{i\Phi/L}$. Concretely, in the limit of +$\Phi/L \to 0$, we have: $$\frac{E[\Phi] - E[0]}{L} \approx \frac{1}{2} \Upsilon(L) \bigg (\frac{\Phi}{L}\bigg)^2 + \cdots$$ -In order to find the ground state under these twisted boundary conditions, we must construct our own variant of the Bose-Hubbard Hamiltonian. Typically you would want to take a peek at the [source code](https://github.com/QuantumKitHub/MPSKitModels.jl/blob/main/src/models/hamiltonians.jl#L379) of `MPSKitModels.jl` to see how these models are defined and tweak it as per your needs. Here we see that applying twisted boundary conditions is equivalent to adding a prefactor of $e^{\pm i\phi}$ in front of the hopping amplitudes. +In order to find the ground state under these twisted boundary conditions, we must construct +our own variant of the Bose-Hubbard Hamiltonian. Typically you would want to take a peek at +the +[source code](https://github.com/QuantumKitHub/MPSKitModels.jl/blob/f4c36d9660a9eab05fa253ffd5c20dc6b7df44cc/src/models/hamiltonians.jl#L379-L409) +of `MPSKitModels.jl` to see how these models are defined and tweak it as per your needs. +Here we see that applying twisted boundary conditions is equivalent to adding a prefactor of +$e^{\pm i\phi}$ in front of the hopping amplitudes. ````julia function bose_hubbard_model_twisted_bc( @@ -323,10 +419,17 @@ superfluid_stiffness_profile(0.01, 0.3, 5, 4) # mott insulator ```` ```@raw html - + ``` -Now that we know what phases to expect, we can plot the phase diagram by scanning over a range of parameters. In general, one could do better by performing a bisection algorithm for each chemical potential to determine the value of the hopping parameter at the transition point, however the 1D Bose-Hubbard model may have two transition points at the same chemical potential which makes this a bit cumbersome to implement robustly. Furthermore, we stick to using the quasi-condensate density as an order parameter since extracting the superfluid density accurately requires a more robust scheme to compute second derivatives which takes us away from the focus of this tutorial. +Now that we know what phases to expect, we can plot the phase diagram by scanning over a +range of parameters. In general, one could do better by performing a bisection algorithm for +each chemical potential to determine the value of the hopping parameter at the transition +point, however the 1D Bose-Hubbard model may have two transition points at the same chemical +potential which makes this a bit cumbersome to implement robustly. Furthermore, we stick to +using the quasi-condensate density as an order parameter since extracting the superfluid +density accurately requires a more robust scheme to compute second derivatives which takes +us away from the focus of this tutorial. ````julia cutoff, D = 4, 10 @@ -347,10 +450,14 @@ heatmap(ts, mus, order_parameters, xlabel = L"t/U", ylabel = L"\mu/U", title = L ```` ```@raw html - + ``` -Although the bond dimension here is quite low, we already see the deformation of the Mott insulator lobes to give way to the well known BKT transition that happens at commensurate density. One can go further and estimate the critical exponents using finite-entanglement scaling procedures on the correlation functions, but these may now be performed with ease using what we have learnt in this tutorial. +Although the bond dimension here is quite low, we already see the deformation of the Mott +insulator lobes to give way to the well known BKT transition that happens at commensurate +density. One can go further and estimate the critical exponents using finite-entanglement +scaling procedures on the correlation functions, but these may now be performed with ease +using what we have learnt in this tutorial. --- diff --git a/docs/src/examples/quantum1d/8.bose-hubbard/main.ipynb b/docs/src/examples/quantum1d/8.bose-hubbard/main.ipynb index e909ec636..13d29b9bd 100644 --- a/docs/src/examples/quantum1d/8.bose-hubbard/main.ipynb +++ b/docs/src/examples/quantum1d/8.bose-hubbard/main.ipynb @@ -19,11 +19,17 @@ "source": [ "# 1D Bose-Hubbard model\n", "\n", - "In this tutorial, we will explore the physics of the one-dimensional Bose–Hubbard model using matrix product states. For the most part, we replicate the results presented in [**Phys. Rev. B 105, 134502**](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.105.134502), which can be consulted for any statements in this tutorial that are not otherwise cited. The Hamiltonian under study is now defined as follows:\n", + "In this tutorial, we will explore the physics of the one-dimensional Bose–Hubbard model\n", + "using matrix product states. For the most part, we replicate the results presented in\n", + "[**Phys. Rev. B 105,\n", + "134502**](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.105.134502), which can be\n", + "consulted for any statements in this tutorial that are not otherwise cited. The Hamiltonian\n", + "under study is now defined as follows:\n", "\n", "$$H = -t \\sum_{i} (\\hat{a}_i^{\\dagger} \\hat{a}_{i+1} + \\hat{a}_{i+1}^{\\dagger} \\hat{a}_i) + \\frac{U}{2} \\sum_i \\hat{n}_i(\\hat{n}_i - 1) - \\mu \\sum_i \\hat{n}_i$$\n", "\n", - "where the bosonic creation and annihilation operators satisfy the canonical commutation relations (CCR):\n", + "where the bosonic creation and annihilation operators satisfy the canonical commutation\n", + "relations (CCR):\n", "\n", "$$[\\hat{a}_i, \\hat{a}_j^{\\dagger}] = \\delta_{ij}.$$\n", "\n", @@ -61,9 +67,27 @@ "\n", "$$\\hat{n} = \\hat{a}^\\dagger \\hat{a} = \\mathrm{diag}(0, 1, 2, \\ldots, n_{\\text{max}}).$$\n", "\n", - "Before moving on, notice that the Hamiltonian is uniform and translationally invariant. Typically, such models are studied on a finite chain of $N$ sites with periodic boundary conditions, but this introduces finite-size effects that are rather annoying to deal with. In contrast, the MPS framework allows us to work directly in the thermodynamic limit, avoiding such artifacts. We will follow this line of exploration in this tutorial and leave finite systems for another time.\n", - "\n", - "In order to work in the thermodynamic limit, we will have to create an `InfiniteMPS`. A complete specification of the MPS requires us to define the physical space and the virtual space of the constituent tensors. At this point is it useful to note that `MPSKit.jl` is powered by `TensorKit.jl` under the hood and has some very generic interfaces in order to allow imposing symmetries of all kinds. As a result, it is sometimes necessary to be a bit more explicit about what we want to do in terms of the vector spaces involved. In this case, we will not consider any symmetries and simply take the most naive approach of working within the `Trivial` sector. The physical space is then `ℂ^(nmax+1)`, and the virtual space is `ℂ^D` where $D$ is some integer chosen to be the bond dimension of the MPS, and `ℂ` is an alias for `ComplexSpace` (typeset as `\\bbC`). As $D$ is increased, one increases the amount of entanglement, i.e, quantum correlations that can be captured by the state." + "Before moving on, notice that the Hamiltonian is uniform and translationally invariant.\n", + "Typically, such models are studied on a finite chain of $N$ sites with periodic boundary\n", + "conditions, but this introduces finite-size effects that are rather annoying to deal with.\n", + "In contrast, the MPS framework allows us to work directly in the thermodynamic limit,\n", + "avoiding such artifacts. We will follow this line of exploration in this tutorial and leave\n", + "finite systems for another time.\n", + "\n", + "In order to work in the thermodynamic limit, we will have to create an\n", + "`InfiniteMPS`. A complete specification of the MPS requires us to define the\n", + "physical space and the virtual space of the constituent tensors. At this point is it useful\n", + "to note that `MPSKit.jl` is powered by\n", + "[`TensorKit.jl`](https://github.com/QuantumKitHub/TensorKit.jl) under the hood and has some\n", + "very generic interfaces in order to allow imposing symmetries of all kinds. As a result, it\n", + "is sometimes necessary to be a bit more explicit about what we want to do in terms of the\n", + "vector spaces involved. In this case, we will not consider any symmetries and simply take\n", + "the most naive approach of working within the `Trivial`\n", + "sector. The physical space is then `ℂ^(nmax+1)`, and the virtual space is `ℂ^D` where $D$ is\n", + "some integer chosen to be the bond dimension of the MPS, and `ℂ` is an alias for\n", + "`ComplexSpace` (typeset as `\\bbC`). As $D$ is increased,\n", + "one increases the amount of entanglement, i.e, quantum correlations that can be captured by\n", + "the state." ], "metadata": {} }, @@ -80,7 +104,12 @@ { "cell_type": "markdown", "source": [ - "This simply initializes a tensor filled with random entries (check out the documentation for other useful constructors). Next, we need the creation and annihilation operators. While we could construct them from scratch, here we will use `MPSKitModels.jl` instead that has predefined operators and models for most well-known lattice models." + "This simply initializes a tensor filled with random entries (check out the documentation for\n", + "other useful constructors). Next, we need the creation and annihilation operators. While we\n", + "could construct them from scratch, here we will use\n", + "[`MPSKitModels.jl`](https://github.com/QuantumKitHub/MPSKitModels.jl) instead that has\n", + "predefined operators and models for most well-known lattice models. In particular, we can\n", + "use `MPSKitModels.a_min` to create the bosonic annihilation operator." ], "metadata": {} }, @@ -98,7 +127,10 @@ { "cell_type": "markdown", "source": [ - "The [] accessor lets us see the underlying array and indeed the operators are exactly what we require. Similarly, the Bose Hubbard model is also predefined (although we will construct our own variant later on)." + "The [] accessor lets us see the underlying array, and indeed the operators are exactly what\n", + "we require. Similarly, the Bose Hubbard model is also predefined in\n", + "`MPSKitModels.bose_hubbard_model` also predefined (although we will construct our\n", + "own variant later on)." ], "metadata": {} }, @@ -114,7 +146,11 @@ { "cell_type": "markdown", "source": [ - "This has created the Hamiltonian operator as a matrix product operator (MPO) which is a convenient form to use in conjunction with MPS. Finally, the ground state optimization may be performed with either `iDMRG` or `VUMPS`. Both should take similar arguments but it is known that VUMPS is typically more efficient for these systems so we proceed with that." + "This has created the Hamiltonian operator as a [matrix product operator](@ref\n", + "InfiniteMPOHamiltonian) (MPO) which is a convenient form to use in conjunction with MPS.\n", + "Finally, the ground state optimization may be performed with either `iDMRG` or\n", + "`VUMPS`. Both should take similar arguments but it is known that VUMPS is typically\n", + "more efficient for these systems so we proceed with that." ], "metadata": {} }, @@ -131,7 +167,9 @@ { "cell_type": "markdown", "source": [ - "This automatically runs the algorithm until a certain error [measure](https://github.com/QuantumKitHub/MPSKit.jl/blob/main/src/algorithms/toolbox.jl#L42) falls below the specified tolerance or the maximum iterations is reached. Let us wrap all this into a convenient function." + "This automatically runs the algorithm until a certain error [measure](@ref\n", + "MPSKit.calc_galerkin) falls below the specified tolerance or the maximum iterations is\n", + "reached. Let us wrap all this into a convenient function." ], "metadata": {} }, @@ -155,7 +193,11 @@ { "cell_type": "markdown", "source": [ - "Now that we have the state, we may compute observables using the `expectation_value` function. It typically expects a `Pair`, `(i1, i2, .., ik) => op` where `op` is a `TensorMap` or `InfiniteMPO` acting over `k` sites. In case of the Hamiltonian, it is not necessary to specify the indices as it spans the whole lattice. We can now plot the correlation function $\\langle \\hat{a}^{\\dagger}_i \\hat{a}_j\\rangle$." + "Now that we have the state, we may compute observables using the `expectation_value`\n", + "function. It typically expects a `Pair`, `(i1, i2, .., ik) => op` where `op` is a\n", + "`TensorMap` or `InfiniteMPO` acting over `k` sites. In case of the Hamiltonian, it is not\n", + "necessary to specify the indices as it spans the whole lattice. We can now plot the\n", + "correlation function $\\langle \\hat{a}^{\\dagger}_i \\hat{a}_j\\rangle$." ], "metadata": {} }, @@ -172,7 +214,8 @@ { "cell_type": "markdown", "source": [ - "We see that the correlations drop off exponentially, indicating the existence of a gapped Mott insulating phase. Let us now shift our parameters to probe other phases." + "We see that the correlations drop off exponentially, indicating the existence of a gapped\n", + "Mott insulating phase. Let us now shift our parameters to probe other phases." ], "metadata": {} }, @@ -191,7 +234,19 @@ { "cell_type": "markdown", "source": [ - "In this case, the correlation function drops off algebraically and eventually saturates as $\\lim_{i \\to \\infty}\\langle\\hat{a}_i^{\\dagger} \\hat{a}_j\\rangle ≈ \\langle \\hat{a}_i^{\\dagger}\\rangle \\langle \\hat{a}_j \\rangle = |\\langle a_i \\rangle|^2 \\neq 0$. This is a signature of long-range order and suggests the existence of a Bose-Einstein condensate. However, this is a bit odd since at zero temperature, the Bose Hubbard model is not expected to break any continuous symmetries ($U(1)$ in this case, corresponding to particle number conservation) due to the Mermin-Wagner theorem. The source of this contradiction lies in the fact that the true 1D superfluid ground state is an extended critical phase exhibiting algebraic decay, however, a finite bond-dimension MPS can only capture exponentially decaying correlations. As a result, the finite bond dimension effectively introduces a length scale into the system in a similar manner as finite-size effects. We can see this clearly by increasing the bond dimension. We also see that the correlation length seems to depend algebraically on the bond dimension as expected from finite-entanglement scaling arguments." + "In this case, the correlation function drops off algebraically and eventually saturates as\n", + "$\\lim_{i \\to \\infty}\\langle\\hat{a}_i^{\\dagger} \\hat{a}_j\\rangle ≈ \\langle \\hat{a}_i^{\\dagger}\\rangle \\langle \\hat{a}_j \\rangle = |\\langle a_i \\rangle|^2 \\neq 0$.\n", + "This is a signature of long-range order and suggests the existence of a Bose-Einstein\n", + "condensate. However, this is a bit odd since at zero temperature, the Bose Hubbard model is\n", + "not expected to break any continuous symmetries ($U(1)$ in this case, corresponding to\n", + "particle number conservation) due to the Mermin-Wagner theorem. The source of this\n", + "contradiction lies in the fact that the true 1D superfluid ground state is an extended\n", + "critical phase exhibiting algebraic decay, however, a finite bond-dimension MPS can only\n", + "capture exponentially decaying correlations. As a result, the finite bond dimension\n", + "effectively introduces a length scale into the system in a similar manner as finite-size\n", + "effects. We can see this clearly by increasing the bond dimension. We also see that the\n", + "correlation length seems to depend algebraically on the bond dimension as expected from\n", + "finite-entanglement scaling arguments." ], "metadata": {} }, @@ -248,7 +303,11 @@ { "cell_type": "markdown", "source": [ - "This shows that any finite bond dimension MPS necessarily breaks the symmetry of the system, forming a Bose-Einstein condensate which introduces erraneous long-distance behaviour of correlation functions. In case of finite bond dimension, it is thus reasonable to associate the finite expectation value of the field operator to the 'quasicondensate' density of the system which vanishes as $D \\to \\infty$." + "This shows that any finite bond dimension MPS necessarily breaks the symmetry of the system,\n", + "forming a Bose-Einstein condensate which introduces erraneous long-distance behaviour of\n", + "correlation functions. In case of finite bond dimension, it is thus reasonable to associate\n", + "the finite expectation value of the field operator to the 'quasicondensate' density of the\n", + "system which vanishes as $D \\to \\infty$." ], "metadata": {} }, @@ -266,14 +325,17 @@ "source": [ "We may now also visualize the momentum distribution function, which is obtained as the Fourier transform of the single-particle density matrix. Starting from the definition of the momentum occupation operator:\n", "\n", - "$$\\hat{a}_k = \\frac{1}{\\sqrt{L}} \\sum_j e^{-ikj} \\hat{a}_j, \\qquad\n", + "$$\n", + "\\hat{a}_k = \\frac{1}{\\sqrt{L}} \\sum_j e^{-ikj} \\hat{a}_j, \\qquad\n", "\\hat{a}_k^\\dagger = \\frac{1}{\\sqrt{L}} \\sum_{j'} e^{ikj'} \\hat{a}_{j'}^\\dagger\n", "$$\n", "\n", "the momentum distribution is\n", "\n", - "$$\\langle \\hat{n}_k \\rangle = \\langle \\hat{a}_k^\\dagger \\hat{a}_k \\rangle\n", - "= \\frac{1}{L} \\sum_{j',j} e^{ik(j'-j)} \\langle \\hat{a}_{j'}^\\dagger \\hat{a}_j \\rangle.$$\n", + "$$\n", + "\\langle \\hat{n}_k \\rangle = \\langle \\hat{a}_k^\\dagger \\hat{a}_k \\rangle\n", + "= \\frac{1}{L} \\sum_{j',j} e^{ik(j'-j)} \\langle \\hat{a}_{j'}^\\dagger \\hat{a}_j \\rangle.\n", + "$$\n", "\n", "For a translationally invariant system, the correlation depends only on the distance $r = j' - j$:\n", "\n", @@ -287,7 +349,11 @@ "\n", "$$\\langle \\hat{n}_k \\rangle = \\sum_{r \\in \\mathbb{Z}} e^{ikr} \\langle \\hat{a}_r^\\dagger \\hat{a}_0 \\rangle$$\n", "\n", - "However, we know that a finite bond dimension MPS introduces a non-zero quasi-condensate density which would give rise to an $\\mathcal{O}(N)$ divergence in the momentum distribution that is not indicative of the true physics of the system. Since we know this contribution vanishes in the infinite bond dimension limit, we instead work with $\\langle \\hat{a}_r^{\\dagger} \\hat{a}_0 \\rangle_c = \\langle \\hat{a}_r^{\\dagger} \\hat{a}_0 \\rangle - |\\langle \\hat{a}\\rangle|^2$." + "However, we know that a finite bond dimension MPS introduces a non-zero quasi-condensate\n", + "density which would give rise to an $\\mathcal{O}(N)$ divergence in the momentum distribution\n", + "that is not indicative of the true physics of the system. Since we know this contribution\n", + "vanishes in the infinite bond dimension limit, we instead work with\n", + "$\\langle \\hat{a}_r^{\\dagger} \\hat{a}_0 \\rangle_c = \\langle \\hat{a}_r^{\\dagger} \\hat{a}_0 \\rangle - |\\langle \\hat{a}\\rangle|^2$." ], "metadata": {} }, @@ -314,13 +380,39 @@ { "cell_type": "markdown", "source": [ - "We see that the density seems to peak around $k=0$, this time seemingly becoming more prominent as $D \\to \\infty$ which seems to suggest again that there is a condensate. However, going by the Penrose-Onsager criterion, the existence of a condensate can be quantified by requiring the leading eigenvalue of the single particle density matrix (i.e, $\\langle \\hat{n}_{k=0}\\rangle = \\sum_j \\langle \\hat{a}_j^{\\dagger} \\hat{a}_0\\rangle$) to diverge as $O(N)$ in the thermodynamic limit. In this case, since the correlations decay as a power law, there is naturally a divergence at low momenta. But this does not imply the existence of a condensate since the order of divergence is much weaker. However, this does indicate the remnants of some kind of condensation in the 1D model despite the quantum fluctuations, leading to the practical utility of defining the concept of a quasicondensate where there is still a notion of phase coherence over short distances.\n", - "\n", - "What this means for us is that, as far as MPS simulations go, we may still utilize the quasicondensate density as an effective order parameter, although it will be less robust as the bond dimension is increased. Alternatively, we realize that the true phase is characterized as being a superfluid (a concept distinct from Bose-Einstein condensation) and can be identified by a non-zero value of the superfluid stiffness (also known as helicity modulus, $\\Upsilon$) as defined by Leggett. Upon applying a phase twist $\\Phi$ to the boundaries of the system, a superfluid phase would suffer an increase in energy whereas an insulating phase would not. In the thermodynamic limit, one could show that the boundary conditions may be considered as periodic and instead uniformly distribute the phase across the chain as $\\hat{a}_i \\to \\hat{a}_i e^{i\\Phi/L}$. Concretely, in the limit of $\\Phi/L \\to 0$, we have:\n", + "We see that the density seems to peak around $k=0$, this time seemingly becoming more\n", + "prominent as $D \\to \\infty$ which seems to suggest again that there is a condensate.\n", + "However, going by the Penrose-Onsager criterion, the existence of a condensate can be\n", + "quantified by requiring the leading eigenvalue of the single particle density matrix (i.e,\n", + "$\\langle \\hat{n}_{k=0}\\rangle = \\sum_j \\langle \\hat{a}_j^{\\dagger} \\hat{a}_0\\rangle$) to\n", + "diverge as $O(N)$ in the thermodynamic limit. In this case, since the correlations decay as\n", + "a power law, there is naturally a divergence at low momenta. But this does not imply the\n", + "existence of a condensate since the order of divergence is much weaker. However, this does\n", + "indicate the remnants of some kind of condensation in the 1D model despite the quantum\n", + "fluctuations, leading to the practical utility of defining the concept of a quasicondensate\n", + "where there is still a notion of phase coherence over short distances.\n", + "\n", + "What this means for us is that, as far as MPS simulations go, we may still utilize the\n", + "quasicondensate density as an effective order parameter, although it will be less robust as\n", + "the bond dimension is increased. Alternatively, we realize that the true phase is\n", + "characterized as being a superfluid (a concept distinct from Bose-Einstein condensation) and\n", + "can be identified by a non-zero value of the superfluid stiffness (also known as helicity\n", + "modulus, $\\Upsilon$) as defined by Leggett. Upon applying a phase twist $\\Phi$ to the\n", + "boundaries of the system, a superfluid phase would suffer an increase in energy whereas an\n", + "insulating phase would not. In the thermodynamic limit, one could show that the boundary\n", + "conditions may be considered as periodic and instead uniformly distribute the phase across\n", + "the chain as $\\hat{a}_i \\to \\hat{a}_i e^{i\\Phi/L}$. Concretely, in the limit of\n", + "$\\Phi/L \\to 0$, we have:\n", "\n", "$$\\frac{E[\\Phi] - E[0]}{L} \\approx \\frac{1}{2} \\Upsilon(L) \\bigg (\\frac{\\Phi}{L}\\bigg)^2 + \\cdots$$\n", "\n", - "In order to find the ground state under these twisted boundary conditions, we must construct our own variant of the Bose-Hubbard Hamiltonian. Typically you would want to take a peek at the [source code](https://github.com/QuantumKitHub/MPSKitModels.jl/blob/main/src/models/hamiltonians.jl#L379) of `MPSKitModels.jl` to see how these models are defined and tweak it as per your needs. Here we see that applying twisted boundary conditions is equivalent to adding a prefactor of $e^{\\pm i\\phi}$ in front of the hopping amplitudes." + "In order to find the ground state under these twisted boundary conditions, we must construct\n", + "our own variant of the Bose-Hubbard Hamiltonian. Typically you would want to take a peek at\n", + "the\n", + "[source code](https://github.com/QuantumKitHub/MPSKitModels.jl/blob/f4c36d9660a9eab05fa253ffd5c20dc6b7df44cc/src/models/hamiltonians.jl#L379-L409)\n", + "of `MPSKitModels.jl` to see how these models are defined and tweak it as per your needs.\n", + "Here we see that applying twisted boundary conditions is equivalent to adding a prefactor of\n", + "$e^{\\pm i\\phi}$ in front of the hopping amplitudes." ], "metadata": {} }, @@ -374,7 +466,14 @@ { "cell_type": "markdown", "source": [ - "Now that we know what phases to expect, we can plot the phase diagram by scanning over a range of parameters. In general, one could do better by performing a bisection algorithm for each chemical potential to determine the value of the hopping parameter at the transition point, however the 1D Bose-Hubbard model may have two transition points at the same chemical potential which makes this a bit cumbersome to implement robustly. Furthermore, we stick to using the quasi-condensate density as an order parameter since extracting the superfluid density accurately requires a more robust scheme to compute second derivatives which takes us away from the focus of this tutorial." + "Now that we know what phases to expect, we can plot the phase diagram by scanning over a\n", + "range of parameters. In general, one could do better by performing a bisection algorithm for\n", + "each chemical potential to determine the value of the hopping parameter at the transition\n", + "point, however the 1D Bose-Hubbard model may have two transition points at the same chemical\n", + "potential which makes this a bit cumbersome to implement robustly. Furthermore, we stick to\n", + "using the quasi-condensate density as an order parameter since extracting the superfluid\n", + "density accurately requires a more robust scheme to compute second derivatives which takes\n", + "us away from the focus of this tutorial." ], "metadata": {} }, @@ -404,7 +503,11 @@ { "cell_type": "markdown", "source": [ - "Although the bond dimension here is quite low, we already see the deformation of the Mott insulator lobes to give way to the well known BKT transition that happens at commensurate density. One can go further and estimate the critical exponents using finite-entanglement scaling procedures on the correlation functions, but these may now be performed with ease using what we have learnt in this tutorial." + "Although the bond dimension here is quite low, we already see the deformation of the Mott\n", + "insulator lobes to give way to the well known BKT transition that happens at commensurate\n", + "density. One can go further and estimate the critical exponents using finite-entanglement\n", + "scaling procedures on the correlation functions, but these may now be performed with ease\n", + "using what we have learnt in this tutorial." ], "metadata": {} }, @@ -424,11 +527,11 @@ "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.11.7" + "version": "1.11.6" }, "kernelspec": { "name": "julia-1.11", - "display_name": "Julia 1.11.7", + "display_name": "Julia 1.11.6", "language": "julia" } }, diff --git a/examples/Cache.toml b/examples/Cache.toml index a4a34c75a..e5f3b32a9 100644 --- a/examples/Cache.toml +++ b/examples/Cache.toml @@ -5,7 +5,7 @@ "2.haldane" = "804433b690faa1ce268a430edb4185af77a38de7a9f53e097795688a53c30c5e" "6.hubbard" = "af14e1df11b2392a69260ce061a716370a2ce50b5c907f0a7d2548e796d543fe" "7.xy-finiteT" = "7afb26fb9ff8ef84722fee3442eaed2ddffda4a5f70a7844b61ce5178093706c" -"8.bose-hubbard" = "63c8df67cede6db14296a87aa86b2e62fde50d0ab5cde6401c07ca692d18ea99" +"8.bose-hubbard" = "665ee9e67ec428b771c0687ccebb232860f7373cfe1ac3081b35032e910dd9fc" "3.ising-dqpt" = "2d314fd05a75c5c91ff81391c7bf166171b35a7b32933e9490756dc584fe8437" "5.haldane-spt" = "3de1b3baa1e4c5dc2252bbe8688d0c6ac8e5d5265deeb6ab66818bbfb02dd4aa" "4.xxz-heisenberg" = "1a2093a182f3ce070b53488484d1024e45496b291c1b74e3f370201d4c056be2" From 9460c2c8ed144ebe1549cbe94dc323b5339a2f8a Mon Sep 17 00:00:00 2001 From: leburgel Date: Wed, 19 Nov 2025 15:11:25 +0100 Subject: [PATCH 6/8] Minor edits --- examples/quantum1d/8.bose-hubbard/main.jl | 40 +++++++++++++---------- 1 file changed, 23 insertions(+), 17 deletions(-) diff --git a/examples/quantum1d/8.bose-hubbard/main.jl b/examples/quantum1d/8.bose-hubbard/main.jl index f525e9ef4..656369727 100644 --- a/examples/quantum1d/8.bose-hubbard/main.jl +++ b/examples/quantum1d/8.bose-hubbard/main.jl @@ -13,7 +13,7 @@ using matrix product states. For the most part, we replicate the results present [**Phys. Rev. B 105, 134502**](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.105.134502), which can be consulted for any statements in this tutorial that are not otherwise cited. The Hamiltonian -under study is now defined as follows: +under study is defined as follows: $$H = -t \sum_{i} (\hat{a}_i^{\dagger} \hat{a}_{i+1} + \hat{a}_{i+1}^{\dagger} \hat{a}_i) + \frac{U}{2} \sum_i \hat{n}_i(\hat{n}_i - 1) - \mu \sum_i \hat{n}_i$$ @@ -61,7 +61,7 @@ Typically, such models are studied on a finite chain of $N$ sites with periodic conditions, but this introduces finite-size effects that are rather annoying to deal with. In contrast, the MPS framework allows us to work directly in the thermodynamic limit, avoiding such artifacts. We will follow this line of exploration in this tutorial and leave -finite systems for another time. +finite systems for another example. In order to work in the thermodynamic limit, we will have to create an [`InfiniteMPS`](@ref). A complete specification of the MPS requires us to define the @@ -99,12 +99,12 @@ md""" The [] accessor lets us see the underlying array, and indeed the operators are exactly what we require. Similarly, the Bose Hubbard model is also predefined in -[`MPSKitModels.bose_hubbard_model`](@extref) also predefined (although we will construct our -own variant later on). +[`MPSKitModels.bose_hubbard_model`](@extref) (although we will construct our own variant +later on). """ -hamiltonian = bose_hubbard_model(InfiniteChain(1), cutoff = cutoff, U = 1, mu = 0.5, t = 0.2) # It is not strictly required to pass InfiniteChain() and is only included for clarity; one may instead pass FiniteChain(N) as well +hamiltonian = bose_hubbard_model(InfiniteChain(1); cutoff = cutoff, U = 1, mu = 0.5, t = 0.2) # It is not strictly required to pass InfiniteChain() and is only included for clarity; one may instead pass FiniteChain(N) as well md""" This has created the Hamiltonian operator as a [matrix product operator](@ref @@ -118,13 +118,13 @@ ground_state, _, _ = find_groundstate(initial_state, hamiltonian, VUMPS(tol = 1. println("Energy: ", expectation_value(ground_state, hamiltonian)) md""" -This automatically runs the algorithm until a certain error [measure](@ref +This automatically runs the algorithm until a certain [error measure](@ref MPSKit.calc_galerkin) falls below the specified tolerance or the maximum iterations is reached. Let us wrap all this into a convenient function. """ function get_ground_state(mu, t, cutoff, D; kwargs...) - hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff = cutoff, U = 1, mu = mu, t = t) + hamiltonian = bose_hubbard_model(InfiniteChain(); cutoff = cutoff, U = 1, mu = mu, t = t) state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) state, _, _ = find_groundstate(state, hamiltonian, VUMPS(; kwargs...)) @@ -161,10 +161,11 @@ $\lim_{i \to \infty}\langle\hat{a}_i^{\dagger} \hat{a}_j\rangle ≈ \langle \hat This is a signature of long-range order and suggests the existence of a Bose-Einstein condensate. However, this is a bit odd since at zero temperature, the Bose Hubbard model is not expected to break any continuous symmetries ($U(1)$ in this case, corresponding to -particle number conservation) due to the Mermin-Wagner theorem. The source of this -contradiction lies in the fact that the true 1D superfluid ground state is an extended -critical phase exhibiting algebraic decay, however, a finite bond-dimension MPS can only -capture exponentially decaying correlations. As a result, the finite bond dimension +particle number conservation) due to the +[Mermin-Wagner theorem](https://en.wikipedia.org/wiki/Mermin%E2%80%93Wagner_theorem). The +source of this contradiction lies in the fact that the true 1D superfluid ground state is an +extended critical phase exhibiting algebraic decay, however, a finite bond-dimension MPS can +only capture exponentially decaying correlations. As a result, the finite bond dimension effectively introduces a length scale into the system in a similar manner as finite-size effects. We can see this clearly by increasing the bond dimension. We also see that the correlation length seems to depend algebraically on the bond dimension as expected from @@ -177,7 +178,7 @@ mu, t = 0.5, 0.2 states = Vector{InfiniteMPS}(undef, length(Ds)) Threads.@threads for idx in eachindex(Ds) - states[idx] = get_ground_state(mu, t, cutoff, Ds[idx], tol = 1.0e-7, verbosity = 1, maxiter = 500) + states[idx] = get_ground_state(mu, t, cutoff, Ds[idx]; tol = 1.0e-7, verbosity = 1, maxiter = 500) end npoints = 400 @@ -225,7 +226,9 @@ system which vanishes as $D \to \infty$. quasicondensate_density = map(state -> abs2(expectation_value(state, (0,) => a_op)), states) md""" -We may now also visualize the momentum distribution function, which is obtained as the Fourier transform of the single-particle density matrix. Starting from the definition of the momentum occupation operator: +We may now also visualize the momentum distribution function, which is obtained as the +Fourier transform of the single-particle density matrix. Starting from the definition of the +momentum occupation operators: ```math \hat{a}_k = \frac{1}{\sqrt{L}} \sum_j e^{-ikj} \hat{a}_j, \qquad @@ -272,7 +275,6 @@ momentum_distribution = vcat(momentum_distribution...)' plot(ks, momentum_distribution, lab = "D = " .* string.(permutedims(Ds)), lw = 1.5, xlabel = "Momentum k", ylabel = L"\langle n_k \rangle", ylim = [0, 50]) md""" - We see that the density seems to peak around $k=0$, this time seemingly becoming more prominent as $D \to \infty$ which seems to suggest again that there is a condensate. However, going by the Penrose-Onsager criterion, the existence of a condensate can be @@ -335,9 +337,13 @@ function superfluid_stiffness_profile(t, mu, D, cutoff, ϵ = 1.0e-4, npoints = 1 energies = zeros(length(phis)) Threads.@threads for idx in eachindex(phis) - hamiltonian_twisted = bose_hubbard_model_twisted_bc(cutoff = cutoff, t = t, mu = mu, U = 1, phi = phis[idx]) + hamiltonian_twisted = bose_hubbard_model_twisted_bc(; + cutoff = cutoff, t = t, mu = mu, U = 1, phi = phis[idx] + ) state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) - state_twisted, _, _ = find_groundstate(state, hamiltonian_twisted, VUMPS(; tol = 1.0e-8, verbosity = 0)) + state_twisted, _, _ = find_groundstate( + state, hamiltonian_twisted, VUMPS(; tol = 1.0e-8, verbosity = 0) + ) energies[idx] = real(expectation_value(state_twisted, hamiltonian_twisted)) end @@ -367,7 +373,7 @@ a_op = a_min(cutoff = cutoff) order_parameters = zeros(length(ts), length(mus)) Threads.@threads for (i, j) in collect(Iterators.product(eachindex(mus), eachindex(ts))) - hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff = cutoff, U = 1, mu = mus[i], t = ts[j]) + hamiltonian = bose_hubbard_model(InfiniteChain(); cutoff = cutoff, U = 1, mu = mus[i], t = ts[j]) init_state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) state, _, _ = find_groundstate(init_state, hamiltonian, VUMPS(; tol = 1.0e-8, verbosity = 0)) order_parameters[i, j] = abs(expectation_value(state, 0 => a_op)) From a80ab3e6c3ff3d367b0c63f0ca6b9310899ec394 Mon Sep 17 00:00:00 2001 From: leburgel Date: Wed, 19 Nov 2025 15:18:30 +0100 Subject: [PATCH 7/8] Break up a few more paragraphs --- examples/quantum1d/8.bose-hubbard/main.jl | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/examples/quantum1d/8.bose-hubbard/main.jl b/examples/quantum1d/8.bose-hubbard/main.jl index 656369727..293a23822 100644 --- a/examples/quantum1d/8.bose-hubbard/main.jl +++ b/examples/quantum1d/8.bose-hubbard/main.jl @@ -22,9 +22,16 @@ relations (CCR): $$[\hat{a}_i, \hat{a}_j^{\dagger}] = \delta_{ij}.$$ -Each lattice site hosts a local Hilbert space corresponding to bosonic occupation states $|n\rangle$, where $(n = 0, 1, 2, \ldots)$. Since this space is formally infinite-dimensional, numerical simulations typically impose a truncation at some maximum occupation number $(n_{\text{max}})$. Such a treatment is justified since it can be observed that the simulation results quickly converge with the cutoff if the filling fraction is kept sufficiently low. +Each lattice site hosts a local Hilbert space corresponding to bosonic occupation states +$|n\rangle$, where $(n = 0, 1, 2, \ldots)$. Since this space is formally +infinite-dimensional, numerical simulations typically impose a truncation at some maximum +occupation number $(n_{\text{max}})$. Such a treatment is justified since it can be observed +that the simulation results quickly converge with the cutoff if the filling fraction is kept +sufficiently low. -Within this truncated space, the local creation and annihilation operators are represented by finite-dimensional matrices. For example, with cutoff $n_{\text{max}}$, the annihilation operator takes the form +Within this truncated space, the local creation and annihilation operators are represented +by finite-dimensional matrices. For example, with cutoff $n_{\text{max}}$, the annihilation +operator takes the form ```math \hat{a} = @@ -242,7 +249,8 @@ the momentum distribution is = \frac{1}{L} \sum_{j',j} e^{ik(j'-j)} \langle \hat{a}_{j'}^\dagger \hat{a}_j \rangle. ``` -For a translationally invariant system, the correlation depends only on the distance $r = j' - j$: +For a translationally invariant system, the correlation depends only on the distance +$r = j' - j$: $$\langle \hat{a}_{j'}^\dagger \hat{a}_j \rangle = C(r) = \langle \hat{a}_r^\dagger \hat{a}_0 \rangle.$$ From f8c9b98e85814290fc9708138da7833ae4d0f773 Mon Sep 17 00:00:00 2001 From: leburgel Date: Wed, 19 Nov 2025 15:51:42 +0100 Subject: [PATCH 8/8] Add renders --- .../quantum1d/8.bose-hubbard/index.md | 83 +++++++++++-------- .../quantum1d/8.bose-hubbard/main.ipynb | 53 +++++++----- examples/Cache.toml | 2 +- 3 files changed, 84 insertions(+), 54 deletions(-) diff --git a/docs/src/examples/quantum1d/8.bose-hubbard/index.md b/docs/src/examples/quantum1d/8.bose-hubbard/index.md index a0ed93336..476461d7a 100644 --- a/docs/src/examples/quantum1d/8.bose-hubbard/index.md +++ b/docs/src/examples/quantum1d/8.bose-hubbard/index.md @@ -22,7 +22,7 @@ using matrix product states. For the most part, we replicate the results present [**Phys. Rev. B 105, 134502**](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.105.134502), which can be consulted for any statements in this tutorial that are not otherwise cited. The Hamiltonian -under study is now defined as follows: +under study is defined as follows: $$H = -t \sum_{i} (\hat{a}_i^{\dagger} \hat{a}_{i+1} + \hat{a}_{i+1}^{\dagger} \hat{a}_i) + \frac{U}{2} \sum_i \hat{n}_i(\hat{n}_i - 1) - \mu \sum_i \hat{n}_i$$ @@ -31,9 +31,16 @@ relations (CCR): $$[\hat{a}_i, \hat{a}_j^{\dagger}] = \delta_{ij}.$$ -Each lattice site hosts a local Hilbert space corresponding to bosonic occupation states $|n\rangle$, where $(n = 0, 1, 2, \ldots)$. Since this space is formally infinite-dimensional, numerical simulations typically impose a truncation at some maximum occupation number $(n_{\text{max}})$. Such a treatment is justified since it can be observed that the simulation results quickly converge with the cutoff if the filling fraction is kept sufficiently low. +Each lattice site hosts a local Hilbert space corresponding to bosonic occupation states +$|n\rangle$, where $(n = 0, 1, 2, \ldots)$. Since this space is formally +infinite-dimensional, numerical simulations typically impose a truncation at some maximum +occupation number $(n_{\text{max}})$. Such a treatment is justified since it can be observed +that the simulation results quickly converge with the cutoff if the filling fraction is kept +sufficiently low. -Within this truncated space, the local creation and annihilation operators are represented by finite-dimensional matrices. For example, with cutoff $n_{\text{max}}$, the annihilation operator takes the form +Within this truncated space, the local creation and annihilation operators are represented +by finite-dimensional matrices. For example, with cutoff $n_{\text{max}}$, the annihilation +operator takes the form ```math \hat{a} = @@ -70,7 +77,7 @@ Typically, such models are studied on a finite chain of $N$ sites with periodic conditions, but this introduces finite-size effects that are rather annoying to deal with. In contrast, the MPS framework allows us to work directly in the thermodynamic limit, avoiding such artifacts. We will follow this line of exploration in this tutorial and leave -finite systems for another time. +finite systems for another example. In order to work in the thermodynamic limit, we will have to create an [`InfiniteMPS`](@ref). A complete specification of the MPS requires us to define the @@ -117,11 +124,11 @@ display((a_op' * a_op)[]) The [] accessor lets us see the underlying array, and indeed the operators are exactly what we require. Similarly, the Bose Hubbard model is also predefined in -[`MPSKitModels.bose_hubbard_model`](@extref) also predefined (although we will construct our -own variant later on). +[`MPSKitModels.bose_hubbard_model`](@extref) (although we will construct our own variant +later on). ````julia -hamiltonian = bose_hubbard_model(InfiniteChain(1), cutoff = cutoff, U = 1, mu = 0.5, t = 0.2) # It is not strictly required to pass InfiniteChain() and is only included for clarity; one may instead pass FiniteChain(N) as well +hamiltonian = bose_hubbard_model(InfiniteChain(1); cutoff = cutoff, U = 1, mu = 0.5, t = 0.2) # It is not strictly required to pass InfiniteChain() and is only included for clarity; one may instead pass FiniteChain(N) as well ```` ```` @@ -146,19 +153,19 @@ println("Energy: ", expectation_value(ground_state, hamiltonian)) ```` ```` -[ Info: VUMPS init: obj = +7.210229722103e-01 err = 5.9367e-01 -[ Info: VUMPS conv 47: obj = -6.757777651163e-01 err = 8.0696299728e-07 time = 0.25 sec -Energy: -0.6757777651163069 + 9.5509045594669e-18im +[ Info: VUMPS init: obj = +5.102229926686e-01 err = 6.1730e-01 +[ Info: VUMPS conv 44: obj = -6.757777651150e-01 err = 9.9483095620e-07 time = 13.31 sec +Energy: -0.6757777651149829 + 1.3734202787398213e-17im ```` -This automatically runs the algorithm until a certain error [measure](@ref +This automatically runs the algorithm until a certain [error measure](@ref MPSKit.calc_galerkin) falls below the specified tolerance or the maximum iterations is reached. Let us wrap all this into a convenient function. ````julia function get_ground_state(mu, t, cutoff, D; kwargs...) - hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff = cutoff, U = 1, mu = mu, t = t) + hamiltonian = bose_hubbard_model(InfiniteChain(); cutoff = cutoff, U = 1, mu = mu, t = t) state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) state, _, _ = find_groundstate(state, hamiltonian, VUMPS(; kwargs...)) @@ -190,7 +197,7 @@ hline!([abs2(expectation_value(ground_state, (0,) => a_op))], ls = :dash, c = :b ```` ```@raw html - + ``` We see that the correlations drop off exponentially, indicating the existence of a gapped @@ -212,10 +219,11 @@ $\lim_{i \to \infty}\langle\hat{a}_i^{\dagger} \hat{a}_j\rangle ≈ \langle \hat This is a signature of long-range order and suggests the existence of a Bose-Einstein condensate. However, this is a bit odd since at zero temperature, the Bose Hubbard model is not expected to break any continuous symmetries ($U(1)$ in this case, corresponding to -particle number conservation) due to the Mermin-Wagner theorem. The source of this -contradiction lies in the fact that the true 1D superfluid ground state is an extended -critical phase exhibiting algebraic decay, however, a finite bond-dimension MPS can only -capture exponentially decaying correlations. As a result, the finite bond dimension +particle number conservation) due to the +[Mermin-Wagner theorem](https://en.wikipedia.org/wiki/Mermin%E2%80%93Wagner_theorem). The +source of this contradiction lies in the fact that the true 1D superfluid ground state is an +extended critical phase exhibiting algebraic decay, however, a finite bond-dimension MPS can +only capture exponentially decaying correlations. As a result, the finite bond dimension effectively introduces a length scale into the system in a similar manner as finite-size effects. We can see this clearly by increasing the bond dimension. We also see that the correlation length seems to depend algebraically on the bond dimension as expected from @@ -228,7 +236,7 @@ mu, t = 0.5, 0.2 states = Vector{InfiniteMPS}(undef, length(Ds)) Threads.@threads for idx in eachindex(Ds) - states[idx] = get_ground_state(mu, t, cutoff, Ds[idx], tol = 1.0e-7, verbosity = 1, maxiter = 500) + states[idx] = get_ground_state(mu, t, cutoff, Ds[idx]; tol = 1.0e-7, verbosity = 1, maxiter = 500) end npoints = 400 @@ -267,7 +275,7 @@ scatter!( ```` ```@raw html - + ``` This shows that any finite bond dimension MPS necessarily breaks the symmetry of the system, @@ -282,16 +290,18 @@ quasicondensate_density = map(state -> abs2(expectation_value(state, (0,) => a_o ```` 7-element Vector{Float64}: - 0.3097427747270374 - 0.28814773216310047 - 0.2701999836479986 - 0.25712727854071404 - 0.24685385462377543 - 0.23539744783500957 - 0.2279965610349862 + 0.30974277868294353 + 0.2881478518741761 + 0.2702001385414079 + 0.2571272814147897 + 0.24685385675266389 + 0.23539778700973993 + 0.2279965507292309 ```` -We may now also visualize the momentum distribution function, which is obtained as the Fourier transform of the single-particle density matrix. Starting from the definition of the momentum occupation operator: +We may now also visualize the momentum distribution function, which is obtained as the +Fourier transform of the single-particle density matrix. Starting from the definition of the +momentum occupation operators: ```math \hat{a}_k = \frac{1}{\sqrt{L}} \sum_j e^{-ikj} \hat{a}_j, \qquad @@ -305,7 +315,8 @@ the momentum distribution is = \frac{1}{L} \sum_{j',j} e^{ik(j'-j)} \langle \hat{a}_{j'}^\dagger \hat{a}_j \rangle. ``` -For a translationally invariant system, the correlation depends only on the distance $r = j' - j$: +For a translationally invariant system, the correlation depends only on the distance +$r = j' - j$: $$\langle \hat{a}_{j'}^\dagger \hat{a}_j \rangle = C(r) = \langle \hat{a}_r^\dagger \hat{a}_0 \rangle.$$ @@ -339,7 +350,7 @@ plot(ks, momentum_distribution, lab = "D = " .* string.(permutedims(Ds)), lw = 1 ```` ```@raw html - + ``` We see that the density seems to peak around $k=0$, this time seemingly becoming more @@ -404,9 +415,13 @@ function superfluid_stiffness_profile(t, mu, D, cutoff, ϵ = 1.0e-4, npoints = 1 energies = zeros(length(phis)) Threads.@threads for idx in eachindex(phis) - hamiltonian_twisted = bose_hubbard_model_twisted_bc(cutoff = cutoff, t = t, mu = mu, U = 1, phi = phis[idx]) + hamiltonian_twisted = bose_hubbard_model_twisted_bc(; + cutoff = cutoff, t = t, mu = mu, U = 1, phi = phis[idx] + ) state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) - state_twisted, _, _ = find_groundstate(state, hamiltonian_twisted, VUMPS(; tol = 1.0e-8, verbosity = 0)) + state_twisted, _, _ = find_groundstate( + state, hamiltonian_twisted, VUMPS(; tol = 1.0e-8, verbosity = 0) + ) energies[idx] = real(expectation_value(state_twisted, hamiltonian_twisted)) end @@ -419,7 +434,7 @@ superfluid_stiffness_profile(0.01, 0.3, 5, 4) # mott insulator ```` ```@raw html - + ``` Now that we know what phases to expect, we can plot the phase diagram by scanning over a @@ -440,7 +455,7 @@ a_op = a_min(cutoff = cutoff) order_parameters = zeros(length(ts), length(mus)) Threads.@threads for (i, j) in collect(Iterators.product(eachindex(mus), eachindex(ts))) - hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff = cutoff, U = 1, mu = mus[i], t = ts[j]) + hamiltonian = bose_hubbard_model(InfiniteChain(); cutoff = cutoff, U = 1, mu = mus[i], t = ts[j]) init_state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D) state, _, _ = find_groundstate(init_state, hamiltonian, VUMPS(; tol = 1.0e-8, verbosity = 0)) order_parameters[i, j] = abs(expectation_value(state, 0 => a_op)) @@ -450,7 +465,7 @@ heatmap(ts, mus, order_parameters, xlabel = L"t/U", ylabel = L"\mu/U", title = L ```` ```@raw html - + ``` Although the bond dimension here is quite low, we already see the deformation of the Mott diff --git a/docs/src/examples/quantum1d/8.bose-hubbard/main.ipynb b/docs/src/examples/quantum1d/8.bose-hubbard/main.ipynb index 13d29b9bd..83681551a 100644 --- a/docs/src/examples/quantum1d/8.bose-hubbard/main.ipynb +++ b/docs/src/examples/quantum1d/8.bose-hubbard/main.ipynb @@ -24,7 +24,7 @@ "[**Phys. Rev. B 105,\n", "134502**](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.105.134502), which can be\n", "consulted for any statements in this tutorial that are not otherwise cited. The Hamiltonian\n", - "under study is now defined as follows:\n", + "under study is defined as follows:\n", "\n", "$$H = -t \\sum_{i} (\\hat{a}_i^{\\dagger} \\hat{a}_{i+1} + \\hat{a}_{i+1}^{\\dagger} \\hat{a}_i) + \\frac{U}{2} \\sum_i \\hat{n}_i(\\hat{n}_i - 1) - \\mu \\sum_i \\hat{n}_i$$\n", "\n", @@ -33,9 +33,16 @@ "\n", "$$[\\hat{a}_i, \\hat{a}_j^{\\dagger}] = \\delta_{ij}.$$\n", "\n", - "Each lattice site hosts a local Hilbert space corresponding to bosonic occupation states $|n\\rangle$, where $(n = 0, 1, 2, \\ldots)$. Since this space is formally infinite-dimensional, numerical simulations typically impose a truncation at some maximum occupation number $(n_{\\text{max}})$. Such a treatment is justified since it can be observed that the simulation results quickly converge with the cutoff if the filling fraction is kept sufficiently low.\n", + "Each lattice site hosts a local Hilbert space corresponding to bosonic occupation states\n", + "$|n\\rangle$, where $(n = 0, 1, 2, \\ldots)$. Since this space is formally\n", + "infinite-dimensional, numerical simulations typically impose a truncation at some maximum\n", + "occupation number $(n_{\\text{max}})$. Such a treatment is justified since it can be observed\n", + "that the simulation results quickly converge with the cutoff if the filling fraction is kept\n", + "sufficiently low.\n", "\n", - "Within this truncated space, the local creation and annihilation operators are represented by finite-dimensional matrices. For example, with cutoff $n_{\\text{max}}$, the annihilation operator takes the form\n", + "Within this truncated space, the local creation and annihilation operators are represented\n", + "by finite-dimensional matrices. For example, with cutoff $n_{\\text{max}}$, the annihilation\n", + "operator takes the form\n", "\n", "$$\n", "\\hat{a} =\n", @@ -72,7 +79,7 @@ "conditions, but this introduces finite-size effects that are rather annoying to deal with.\n", "In contrast, the MPS framework allows us to work directly in the thermodynamic limit,\n", "avoiding such artifacts. We will follow this line of exploration in this tutorial and leave\n", - "finite systems for another time.\n", + "finite systems for another example.\n", "\n", "In order to work in the thermodynamic limit, we will have to create an\n", "`InfiniteMPS`. A complete specification of the MPS requires us to define the\n", @@ -129,8 +136,8 @@ "source": [ "The [] accessor lets us see the underlying array, and indeed the operators are exactly what\n", "we require. Similarly, the Bose Hubbard model is also predefined in\n", - "`MPSKitModels.bose_hubbard_model` also predefined (although we will construct our\n", - "own variant later on)." + "`MPSKitModels.bose_hubbard_model` (although we will construct our own variant\n", + "later on)." ], "metadata": {} }, @@ -138,7 +145,7 @@ "outputs": [], "cell_type": "code", "source": [ - "hamiltonian = bose_hubbard_model(InfiniteChain(1), cutoff = cutoff, U = 1, mu = 0.5, t = 0.2) # It is not strictly required to pass InfiniteChain() and is only included for clarity; one may instead pass FiniteChain(N) as well" + "hamiltonian = bose_hubbard_model(InfiniteChain(1); cutoff = cutoff, U = 1, mu = 0.5, t = 0.2) # It is not strictly required to pass InfiniteChain() and is only included for clarity; one may instead pass FiniteChain(N) as well" ], "metadata": {}, "execution_count": null @@ -167,7 +174,7 @@ { "cell_type": "markdown", "source": [ - "This automatically runs the algorithm until a certain error [measure](@ref\n", + "This automatically runs the algorithm until a certain [error measure](@ref\n", "MPSKit.calc_galerkin) falls below the specified tolerance or the maximum iterations is\n", "reached. Let us wrap all this into a convenient function." ], @@ -178,7 +185,7 @@ "cell_type": "code", "source": [ "function get_ground_state(mu, t, cutoff, D; kwargs...)\n", - " hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff = cutoff, U = 1, mu = mu, t = t)\n", + " hamiltonian = bose_hubbard_model(InfiniteChain(); cutoff = cutoff, U = 1, mu = mu, t = t)\n", " state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D)\n", " state, _, _ = find_groundstate(state, hamiltonian, VUMPS(; kwargs...))\n", "\n", @@ -239,10 +246,11 @@ "This is a signature of long-range order and suggests the existence of a Bose-Einstein\n", "condensate. However, this is a bit odd since at zero temperature, the Bose Hubbard model is\n", "not expected to break any continuous symmetries ($U(1)$ in this case, corresponding to\n", - "particle number conservation) due to the Mermin-Wagner theorem. The source of this\n", - "contradiction lies in the fact that the true 1D superfluid ground state is an extended\n", - "critical phase exhibiting algebraic decay, however, a finite bond-dimension MPS can only\n", - "capture exponentially decaying correlations. As a result, the finite bond dimension\n", + "particle number conservation) due to the\n", + "[Mermin-Wagner theorem](https://en.wikipedia.org/wiki/Mermin%E2%80%93Wagner_theorem). The\n", + "source of this contradiction lies in the fact that the true 1D superfluid ground state is an\n", + "extended critical phase exhibiting algebraic decay, however, a finite bond-dimension MPS can\n", + "only capture exponentially decaying correlations. As a result, the finite bond dimension\n", "effectively introduces a length scale into the system in a similar manner as finite-size\n", "effects. We can see this clearly by increasing the bond dimension. We also see that the\n", "correlation length seems to depend algebraically on the bond dimension as expected from\n", @@ -260,7 +268,7 @@ "states = Vector{InfiniteMPS}(undef, length(Ds))\n", "\n", "Threads.@threads for idx in eachindex(Ds)\n", - " states[idx] = get_ground_state(mu, t, cutoff, Ds[idx], tol = 1.0e-7, verbosity = 1, maxiter = 500)\n", + " states[idx] = get_ground_state(mu, t, cutoff, Ds[idx]; tol = 1.0e-7, verbosity = 1, maxiter = 500)\n", "end\n", "\n", "npoints = 400\n", @@ -323,7 +331,9 @@ { "cell_type": "markdown", "source": [ - "We may now also visualize the momentum distribution function, which is obtained as the Fourier transform of the single-particle density matrix. Starting from the definition of the momentum occupation operator:\n", + "We may now also visualize the momentum distribution function, which is obtained as the\n", + "Fourier transform of the single-particle density matrix. Starting from the definition of the\n", + "momentum occupation operators:\n", "\n", "$$\n", "\\hat{a}_k = \\frac{1}{\\sqrt{L}} \\sum_j e^{-ikj} \\hat{a}_j, \\qquad\n", @@ -337,7 +347,8 @@ "= \\frac{1}{L} \\sum_{j',j} e^{ik(j'-j)} \\langle \\hat{a}_{j'}^\\dagger \\hat{a}_j \\rangle.\n", "$$\n", "\n", - "For a translationally invariant system, the correlation depends only on the distance $r = j' - j$:\n", + "For a translationally invariant system, the correlation depends only on the distance\n", + "$r = j' - j$:\n", "\n", "$$\\langle \\hat{a}_{j'}^\\dagger \\hat{a}_j \\rangle = C(r) = \\langle \\hat{a}_r^\\dagger \\hat{a}_0 \\rangle.$$\n", "\n", @@ -447,9 +458,13 @@ " energies = zeros(length(phis))\n", "\n", " Threads.@threads for idx in eachindex(phis)\n", - " hamiltonian_twisted = bose_hubbard_model_twisted_bc(cutoff = cutoff, t = t, mu = mu, U = 1, phi = phis[idx])\n", + " hamiltonian_twisted = bose_hubbard_model_twisted_bc(;\n", + " cutoff = cutoff, t = t, mu = mu, U = 1, phi = phis[idx]\n", + " )\n", " state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D)\n", - " state_twisted, _, _ = find_groundstate(state, hamiltonian_twisted, VUMPS(; tol = 1.0e-8, verbosity = 0))\n", + " state_twisted, _, _ = find_groundstate(\n", + " state, hamiltonian_twisted, VUMPS(; tol = 1.0e-8, verbosity = 0)\n", + " )\n", " energies[idx] = real(expectation_value(state_twisted, hamiltonian_twisted))\n", " end\n", "\n", @@ -489,7 +504,7 @@ "order_parameters = zeros(length(ts), length(mus))\n", "\n", "Threads.@threads for (i, j) in collect(Iterators.product(eachindex(mus), eachindex(ts)))\n", - " hamiltonian = bose_hubbard_model(InfiniteChain(), cutoff = cutoff, U = 1, mu = mus[i], t = ts[j])\n", + " hamiltonian = bose_hubbard_model(InfiniteChain(); cutoff = cutoff, U = 1, mu = mus[i], t = ts[j])\n", " init_state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D)\n", " state, _, _ = find_groundstate(init_state, hamiltonian, VUMPS(; tol = 1.0e-8, verbosity = 0))\n", " order_parameters[i, j] = abs(expectation_value(state, 0 => a_op))\n", diff --git a/examples/Cache.toml b/examples/Cache.toml index e5f3b32a9..0c2fdb60b 100644 --- a/examples/Cache.toml +++ b/examples/Cache.toml @@ -5,7 +5,7 @@ "2.haldane" = "804433b690faa1ce268a430edb4185af77a38de7a9f53e097795688a53c30c5e" "6.hubbard" = "af14e1df11b2392a69260ce061a716370a2ce50b5c907f0a7d2548e796d543fe" "7.xy-finiteT" = "7afb26fb9ff8ef84722fee3442eaed2ddffda4a5f70a7844b61ce5178093706c" -"8.bose-hubbard" = "665ee9e67ec428b771c0687ccebb232860f7373cfe1ac3081b35032e910dd9fc" +"8.bose-hubbard" = "4ca414d4b76bbffbd88a7c153dce245c9fab5b7cee0fe0b706b9a10ed87e29e9" "3.ising-dqpt" = "2d314fd05a75c5c91ff81391c7bf166171b35a7b32933e9490756dc584fe8437" "5.haldane-spt" = "3de1b3baa1e4c5dc2252bbe8688d0c6ac8e5d5265deeb6ab66818bbfb02dd4aa" "4.xxz-heisenberg" = "1a2093a182f3ce070b53488484d1024e45496b291c1b74e3f370201d4c056be2"