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/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..476461d7a
--- /dev/null
+++ b/docs/src/examples/quantum1d/8.bose-hubbard/index.md
@@ -0,0 +1,480 @@
+```@meta
+EditURL = "../../../../../examples/quantum1d/8.bose-hubbard/main.jl"
+```
+
+[](https://mybinder.org/v2/gh/QuantumKitHub/MPSKit.jl/gh-pages?filepath=dev/examples/quantum1d/8.bose-hubbard/main.ipynb)
+[](https://nbviewer.jupyter.org/github/QuantumKitHub/MPSKit.jl/blob/gh-pages/dev/examples/quantum1d/8.bose-hubbard/main.ipynb)
+[](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 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
+
+```math
+\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,
+
+```math
+\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 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
+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
+initial_state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D)
+````
+
+````
+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`](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
+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 in
+[`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
+````
+
+````
+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](@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))
+println("Energy: ", expectation_value(ground_state, hamiltonian))
+````
+
+````
+[ 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
+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)
+ 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 = 1.0e-6, verbosity = 2, maxiter = 500)
+````
+
+````
+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`](@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)
+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 = 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)
+````
+
+```@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](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
+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 = 1.0e-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.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 operators:
+
+```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
+
+```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$:
+
+$$\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/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(
+ 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)))
+
+ 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})
+ 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, ϵ = 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]
+ )
+ state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D)
+ state_twisted, _, _ = find_groundstate(
+ state, hamiltonian_twisted, VUMPS(; tol = 1.0e-8, verbosity = 0)
+ )
+ energies[idx] = real(expectation_value(state_twisted, hamiltonian_twisted))
+ end
+
+ 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
+
+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 = 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")
+````
+
+```@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..83681551a
--- /dev/null
+++ b/docs/src/examples/quantum1d/8.bose-hubbard/main.ipynb
@@ -0,0 +1,554 @@
+{
+ "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\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 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\n",
+ "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",
+ "$|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\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",
+ "\\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.\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 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",
+ "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": {}
+ },
+ {
+ "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\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": {}
+ },
+ {
+ "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\n",
+ "we require. Similarly, the Bose Hubbard model is also predefined in\n",
+ "`MPSKitModels.bose_hubbard_model` (although we will construct our own variant\n",
+ "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](@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": {}
+ },
+ {
+ "outputs": [],
+ "cell_type": "code",
+ "source": [
+ "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": {},
+ "execution_count": null
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "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": {}
+ },
+ {
+ "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 = 1.0e-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`\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": {}
+ },
+ {
+ "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\n",
+ "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 = 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)"
+ ],
+ "metadata": {},
+ "execution_count": null
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "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\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",
+ "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 = 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",
+ "\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(\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
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "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": {}
+ },
+ {
+ "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\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",
+ "\\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",
+ "$$\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\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",
+ "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\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": {}
+ },
+ {
+ "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\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])"
+ ],
+ "metadata": {},
+ "execution_count": null
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "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\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": {}
+ },
+ {
+ "outputs": [],
+ "cell_type": "code",
+ "source": [
+ "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",
+ "\n",
+ " interaction_term = N * (N - id(domain(N)))\n",
+ "\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",
+ " 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, ϵ = 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(;\n",
+ " cutoff = cutoff, t = t, mu = mu, U = 1, phi = phis[idx]\n",
+ " )\n",
+ " state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D)\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",
+ " 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",
+ "\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\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": {}
+ },
+ {
+ "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 = 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\")"
+ ],
+ "metadata": {},
+ "execution_count": null
+ },
+ {
+ "cell_type": "markdown",
+ "source": [
+ "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": {}
+ },
+ {
+ "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.6"
+ },
+ "kernelspec": {
+ "name": "julia-1.11",
+ "display_name": "Julia 1.11.6",
+ "language": "julia"
+ }
+ },
+ "nbformat": 4
+}
\ No newline at end of file
diff --git a/examples/Cache.toml b/examples/Cache.toml
index 8820b4277..0c2fdb60b 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" = "4ca414d4b76bbffbd88a7c153dce245c9fab5b7cee0fe0b706b9a10ed87e29e9"
"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..293a23822
--- /dev/null
+++ b/examples/quantum1d/8.bose-hubbard/main.jl
@@ -0,0 +1,398 @@
+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 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
+
+```math
+\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,
+
+```math
+\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 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
+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`](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
+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 in
+[`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
+
+md"""
+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](@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)
+ 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 = 1.0e-6, verbosity = 2, maxiter = 500)
+
+md"""
+
+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.
+"""
+
+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)
+
+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](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
+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 = 1.0e-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 operators:
+
+```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
+
+```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$:
+
+$$\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/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(
+ 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)))
+
+ 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})
+ 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, ϵ = 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]
+ )
+ state = InfiniteMPS(ℂ^(cutoff + 1), ℂ^D)
+ state_twisted, _, _ = find_groundstate(
+ state, hamiltonian_twisted, VUMPS(; tol = 1.0e-8, verbosity = 0)
+ )
+ energies[idx] = real(expectation_value(state_twisted, hamiltonian_twisted))
+ end
+
+ 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
+
+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 = 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")
+
+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.
+"""
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(