GitHub

Conventions and the TensorMap type

This package uses of the TensorMap data structure exported by the TensorKit package. You should read the relevant TensorKit documentation for details on the implementation of this data structure, however here we provide a brief summary of the essentials relevant to this packages function. When using TimeEvolutionPEPO, you should always load TensorKit alongside:

using TimeEvolutionPEPO, TensorKit

Conventions

TimeEvolutionPEPO follows the convention established by TensorRenormalizationGroups that virtual indices are assigned to the domain and physical indices the codomain. Specifically, the ordering of PEPO tensors in this package can be summarized by the following code;

julia> ket = bra = ℂ^2; # Typical phyiscal Hilbert space dimension
julia> east = south = west = north = ℂ^4; # Example virtual dimension
julia> pepo = rand(ket * bra', east * south * west' * north');
julia> space(pepo)(ℂ^2 ⊗ (ℂ^2)') ← (ℂ^4 ⊗ ℂ^4 ⊗ (ℂ^4)' ⊗ (ℂ^4)')

Note the dual spaces on the bra, west, and north vector spaces.

Warn

The convention here is to order the ket space before the bra space such that

\[A \rho B \to (A \otimes B^{\top}) | \rho \rangle\rangle\]

where $| \rho \rangle\rangle$ is the vector obtained by concatenating the rows together. Julia, being a language with column-major ordering, employs the opposite convention with respect to the reshape function:

julia> ρ = ["a" "b"; "c" "d"]2×2 Matrix{String}:
 "a"  "b"
 "c"  "d"
julia> reshape(ρ, 4)4-element Vector{String}:
 "a"
 "c"
 "b"
 "d"

In terms of local operator transformations:

julia> ρ = rand(2,2); A = rand(2,2); B = rand(2,2);
julia> reshape(kron(B', A) * reshape(ρ, 4), 2, 2) ≈ A * ρ * Btrue
julia> reshape(kron(A, B') * reshape(ρ, 4), 2, 2) ≈ A * ρ * Bfalse

This is worth bearing in mind if you wish to import a vectorized density matrix obtained outside of this package.

Conversion to and from a dense array

If you wish to convert a TensorMap to an Array, then it is straightforward to do so:

julia> s = ℂ^2ℂ^2
julia> convert(Array, TensorMap(rand, s * s, s))2×2×2 Array{Float64, 3}:
[:, :, 1] =
 0.849073  0.965757
 0.550203  0.520191

[:, :, 2] =
 0.0155905  0.0496271
 0.211519   0.695836

however this is not a free operation. Conversely, a TensorMap can be constructed directly by supplying the data as an Array whose size matches the space supplied to the TensorMap constructor:

julia> data = rand(2,2,2)2×2×2 Array{Float64, 3}:
[:, :, 1] =
 0.939664  0.42671
 0.242389  0.4218

[:, :, 2] =
 0.637818  0.208392
 0.133984  0.779115
julia> TensorMap(data, s * s, s)TensorMap((ℂ^2 ⊗ ℂ^2) ← ℂ^2):
[:, :, 1] =
 0.9396639976154892  0.4267097351767837
 0.2423894925462814  0.42180000961623243

[:, :, 2] =
 0.6378175095936836   0.20839191562868242
 0.13398402503890183  0.7791148493682717
julia> TensorMap(data, s, s * s)TensorMap(ℂ^2 ← (ℂ^2 ⊗ ℂ^2)):
[:, :, 1] =
 0.9396639976154892  0.4267097351767837
 0.2423894925462814  0.42180000961623243

[:, :, 2] =
 0.6378175095936836   0.20839191562868242
 0.13398402503890183  0.7791148493682717

Irrep ordering in symmetric spaces

You may wish to make use of a symmetric local vector space for simulations that, for example, preserve total spin or particle number. When constructing symmetric vector spaces, TensorKit prescribes and ordering of the basis vectors that may not align with what you expect. For example, consider a spin-1 particle. You may wish to make use of a possible U(1) symmetry by using writing the basis in terms of the sectors of well-defined $z$-axis spin projection. The typical ordering of the basis is then $\{| 1 \rangle, | 0 \rangle, | -1 \rangle \}$ such that the spin-$z$ operator looks like:

\[S_z = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & -1 \end{pmatrix}\]

However the order assumed by TensorKit is different:

julia> U1Space(1 => 1, 0 => 1, -1 => 1)Rep[TensorKitSectors.U₁](0=>1, 1=>1, -1=>1)

The basis of charge sectors is ordered instead like $\{| 0 \rangle, | 1 \rangle, | -1 \rangle \}$. When constructing a TensorMap with symmetric vector spaces, one should make sure the data adheres to correct charge sector ordering. TimeEvolutionPEPO exports the permutebasis function:

julia> perm = (2,1,3); # Need to swap first two basis vectors in Sz

julia> Sz = SPIN1[3]
3×3 Matrix{ComplexF64}:
 1.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.0+0.0im

julia> permutebasis(Sz, perm)
3×3 Matrix{ComplexF64}:
 0.0+0.0im  0.0+0.0im   0.0+0.0im
 0.0+0.0im  1.0+0.0im   0.0+0.0im
 0.0+0.0im  0.0+0.0im  -1.0+0.0im

to perform the correct change-of-basis transformation.