States and Operators
In the previous chapter we have seen how to define a Hamiltonian for a lattice model. In this chapter we will see how to work with other types of operators, such as measurements and diagonalizing.
Builtins
There are some basic building blocks that can be used to create states and operators. These functions are taken from QuantumOptics.jl package, but extended to work with lattices:
basisstate
creates a state vector in the basis of a lattice.diagonaloperator
creates a diagonal operator in the basis of a lattice.transition
creates a transition operator between two states in the basis of a lattice.number
,destroy
, andcreate
are the number, annihilation, and creation operators in the basis of a lattice.
Remember the general rule of thumb for working with lattices: You can provide a lattice (and optionally the internal basis) instead of a QuantumOptics.Basis
, and a site (or a (site, internal_index)
tuple) instead of an integer index.
For example, the basisstate
function creates a state vector in the basis of a lattice. Normally it takes a QuantumOptics.Basis
and an index as arguments, and returns a Ket
object. However, you can use a lattice and a site as arguments to create a state vector in the basis of the lattice:
julia> using LatticeModels
julia> l = SquareLattice(6, 6);
julia> site = l[!, x = 3, y = 3]
2-dim Bravais lattice site in 2D space at [3.0, 3.0]
julia> l_bas, ind = LatticeBasis(l), site_index(l, site);
julia> psi1 = basisstate(l_bas, ind); # The pure QuantumOptics.jl way
julia> psi2 = basisstate(l, site); # The LatticeModels.jl way
julia> psi1 == psi2
true
If you work in a composite system (e. g. a lattice with a spin degree of freedom), you can either use composite indexing or construct the state as tensor product of the states for each subsystem:
julia> spin = SpinBasis(1//2)
Spin(1/2)
julia> psi1_composite = basisstate(l, spin, (site, 1));
julia> psi2_composite = basisstate(spin, 1) ⊗ basisstate(l, site); # Note the order!
julia> psi1_composite == psi2_composite
true
The order of the tensor product is important. The first argument is the on-site degrees of freedom, and the second is the lattice. This convention is consistent in the rest of the package — the reason behind this is performance of construct_operator
and OperatorBuilder
.
If this order is not followed, you will probably get an error somewhere in your calculations.
The diagonaloperator
function creates a diagonal operator in the basis of a lattice. Normally it takes a QuantumOptics.Basis
and a vector of values (or a single value) as arguments, and returns an Operator
object. There are several convenient ways to use this function with lattices. As an example let's consider the position operator in the basis of a lattice:
julia> xval = coordvalue(l, :x)
LatticeValue{Float64} on a 36-site SquareLattice in 2D space Values stored in a Vector{Float64}: [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 2.0, 2.0, 2.0, 2.0 … 5.0, 5.0, 5.0, 5.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0]
julia> X_1 = diagonaloperator(l_bas, xval.values); # The pure QuantumOptics.jl way
julia> X_2 = diagonaloperator(l, :x); # The LatticeModels.jl way, with a site property
julia> X_3 = diagonaloperator(xval); # The LatticeModels.jl way, with a LatticeValue
julia> X_1 == X_2 == X_3
true
This notation allows converting any LatticeValue
or Site parameter to an operator. Hence, diagonaloperator(l, Coord(1))
is also valid and will return the same operator.
Generally, a good way to create a custom diagonal operator is by using the LatticeValue
approach.
Consider this example: in the Haldane model the diagonal part is $m$ on the A sublattice and $-m$ on the B sublattice. You can create this operator with the following code:
l = HoneycombLattice(6, 6)
m = 3
ms = LatticeValue(l) do site
site.index == 1 ? m : -m
end
Op = diagonaloperator(ms)
Or, the same thing as a one-liner:
l = HoneycombLattice(6, 6)
m = 3
Op = diagonaloperator((LatticeValue(l, :index) .- 1.5) .* 2 .* m)
Also note the coordoperator
and coordoperators
functions that do the same thing as coordvalue
and coordvalues
, but return the operator instead of the value:
julia> X, Y = coordoperators(l);
julia> X == X_1
true
Measurements
The most common type of measurements is the local density: the average number of particles at each site. This can be calculated using the localdensity
function — it takes a state (a QuantumOptics.Ket
vector or a QuantumOptics.Operator
representing the density matrix) and returns a LatticeValue
with the density at each site.
using LatticeModels, Plots
l = GrapheneRibbon(6, 4)
H = tightbinding_hamiltonian(l)
d = localdensity(groundstate(H))
plot(d)
The localdensity
function uses the following formula to calculate the density at each site: $\rho_i = \text{Tr}(\hat{n}_i \hat{\rho})$, where $\hat{n}_i$ is the number operator at site $i$ and $\hat{\rho}$ is the density matrix. Note that if the values are complex, the function will return the real part of them. This is what makes the next example work.
The local Chern marker is a quantity that can be used to detect topological phases in a lattice model. It can be calculated using the following formula:
\[\mathcal{C}(r) = 4\pi \text{Im} \langle r | P X P Y P | r \rangle\]
Here $P$ is the projector onto the occupied states (e. g. the density matrix), and $X$ and $Y$ are the position operators. The local Chern marker is a real number that can be calculated for each site in the lattice.
Let's do this for a QWZ model Hamiltonian on a square lattice:
using LatticeModels, Plots
l = SquareLattice(6, 6)
sys = l ⊗ SpinBasis(1//2)
ms = ones(l)
ms[x = 3 .. 4, y = 3 .. 4] .= -1
H = qwz(l, ms)
P = densitymatrix(H, mu=0, statistics=FermiDirac)
X, Y = coordoperators(sys)
c = localdensity(-4π * im * P * X * P * Y * P)
heatmap(c, title="Local Chern marker")
A generalization of the density measurement is the localexpect
function. It takes a local operator $\hat{A}$ and a state, and calculates the expectation value of the operator at each site: $A_i = \text{Tr}((\hat{A} \otimes \hat{n}_i) \cdot \hat{\rho})$, where $\hat{n}_i$ is the number operator at site $i$ and $\hat{\rho}$ is the density matrix. Note that the result is a complex numbered LatticeValue
.
As an example, let us visualize a spin wavefunction on a square lattice:
using LatticeModels, Plots
l = SquareLattice(10, 10)
x, y = coordvalues(l)
spin = SpinBasis(1//2)
gauss = @. exp(-0.05 * ((x - 5.5) ^ 2 + (y - 5.5) ^ 2))
wave = @. exp(im * (x + y))
ψ = basisstate(spin, 1) ⊗ (@. gauss .* wave) + basisstate(spin, 2) ⊗ (@. gauss * conj(wave))
normalize!(ψ)
σx = sigmax(spin)
p = plot(layout=2, size=(800, 400))
plot!(p[1], localdensity(ψ), title="Local density")
plot!(p[2], localexpect(σx, ψ) .|> real, title="σx projection")
Diagonalizing
To diagonalize a Hamiltonian or any other operator, you can use the diagonalize
function. It takes an operator and returns a EigenSystem
object with the eigenvalues and eigenvectors of the operator.
using LatticeModels, Plots
l = GrapheneRibbon(6, 4)
H = haldane(l, 0.1, 1)
eig = diagonalize(H)
Diagonalized Hamiltonian (48 eigenvectors)
Eigenvalues in range -4.27959 .. 4.27959
System: One particle on 48-site HoneycombLattice in 2D space
This struct simplifies the access to the eigenvalues and eigenvectors of the operator. You can access the eigenvalues with eig.values
, and eigenvectors as Ket
s can be obtained with the bracket notation eig[i]
or eig[value = E]
:
# The states are sorted by real part of the eigenvalues, so
psi = eig[1] # `psi` is the ground state
psi2 = eig[value = 0] # `psi2` is the state with zero energy
p = plot(layout = @layout[a b; c], size=(800, 800))
plot!(p[1], localdensity(psi), title="Ground state")
plot!(p[2], localdensity(psi2), title="Zero energy state")
scatter!(p[3], eig.values, title="Spectrum", lab="")
You can find the ground state of a Hamiltonian in one line using the groundstate
function:
psi = groundstate(H)
To evaluate both the ground state and its energy, use the findgroundstate
function:
E0, psi = findgroundstate(H)
The diagonalize
function under its hood uses the eigen
function from the LinearAlgebra
standard library. However, this does not work for non-trivial matrix types (e. g. sparse matrices or GPU arrays). For such cases you can pass the second argument to the diagonalize
function, which is a Symbol
indicating the method to use for diagonalization. To use the eigsolve
function from the KrylovKit.jl
package, you can pass :krylovkit
as the second argument. Since this solves the eigenvalue problem iteratively, you can also pass the keyword arguments: n
for the number of eigenvalues to compute, v0
for the initial guess, and the keyword arguments for the eigsolve
function.
l = SquareLattice(100, 100) # A really big lattice
V = LatticeValue(l) do site # A potential well
x, y = site
return (x - 50)^2 + (y - 50)^2 < 40^2 ? -0.5 : 0.0
end
H = tightbinding_hamiltonian(l, V)
eig = diagonalize(H, :krylovkit, n=16) # Compute only 9 eigenvalues with smallest real part
p = plot(layout=16, leg=false, size=(1000, 900))
for i in 1:16
histogram2d!(p[i], localdensity(eig[i]), title="E = $(round(eig.values[i], digits=5))")
end
plot!()
EigenSystem
objects have a wide range of applications in this package. One of them is creating equilibrium states for a given Hamiltonian. This can be done using the densitymatrix
function, which is described in the next section.
Density matrix
After you diagonalize a Hamiltonian, you can calculate the density matrix for the system. Use the densitymatrix
function to do this:
using LatticeModels
l = SquareLattice(6, 6)
sys = System(l, SpinBasis(1//2), mu=0, statistics=FermiDirac, T=0.1)
H = tightbinding_hamiltonian(sys)
eig = diagonalize(H)
P1 = densitymatrix(eig) # Use the default parameters from the `System`
P2 = densitymatrix(eig,
statistics=BoseEinstein, T=0, mu=1) # Or you can override them
┌ Info: Creating density matrix: FermiDirac distribution, T = 0.1, μ = 0.0
└ set `info=false` to disable this message
┌ Info: Creating density matrix: BoseEinstein distribution, T = 0, μ = 1
└ set `info=false` to disable this message
Note that the densitymatrix
function can also be applied to a Hamiltonian
object, in which case it will first diagonalize the Hamiltonian and then calculate the density matrix:
P1_1 = densitymatrix(H)
@assert P1 ≈ P1_1
P2_1 = densitymatrix(H, statistics=BoseEinstein, T=0, mu=1)
@assert P2 ≈ P2_1
┌ Info: Creating density matrix: FermiDirac distribution, T = 0.1, μ = 0.0
└ set `info=false` to disable this message
┌ Info: Creating density matrix: BoseEinstein distribution, T = 0, μ = 1
└ set `info=false` to disable this message
The densitymatrix
function uses a simple formula to calculate the density matrix:
\[\hat{\rho} = \sum_{i} \rho(E_i) | \psi_i \rangle \langle \psi_i |\]
Here $E_i$ are the eigenvalues of the operator, $\psi_i$ are the corresponding eigenvectors, and $\rho(E)$ is the distribution function defined by the statistics
, T
, and mu
parameters. By default, when no additional parameters are passed to the System
or densitymatrix
, the density matrix will be a thermal state at zero temperature.
The basis for these computations is the projector
function, which takes a function p
and an EigenSystem
object d
that represents the diagonalized operator $\hat{O}$. The return value is an operator $\hat{P}$ defined by the formula:
\[\hat{O} = \sum_{i} E_i | \psi_i \rangle \langle \psi_i |, \hspace{1cm} \hat{P} = \sum_{i} p(E_i) | \psi_i \rangle \langle \psi_i |\]
Here $E_i$ are the eigenvalues of the operator $\hat{O}$, and $| \psi_i \rangle$ are the corresponding eigenvectors. The function p
is applied to the eigenvalues to obtain the diagonal elements of the density matrix. Here is an example of how to use this function:
l = SquareLattice(6, 6)
H = tightbinding_hamiltonian(l)
eig = diagonalize(H)
P1 = projector(x -> x < 0, eig) # Projector onto the states with energy < 0
P2 = projector(x -> 1 / (1 + exp(x)), eig) # Fermi-Dirac distribution
P3 = projector(eig[1:4]) # Projector onto the first 4 states
# Note how we slice the `eig` object to get the first 4 states