Constructing models

One of the strengths of tensor network methods is that, by working directly with state vectors and matrix representations of time-evolution operators, they do not require deriving specific equations for different models. For this reason, TimeEvolutionPEPO aims to have very flexible support for implementing different models.

Terms, Hamiltonians, and Liouvillians

Models are generated by appending objects that are subtypes of AbstractTerm to an instance of the Model type. Each element in the iterable Model represents a layer in a Suzuki-Trotter decomposition, and is represented by the TrotterLayer type. Multiple terms can appear in the same TrotterLayer indicating that these terms are not trotterised over, i.e. the sum of these terms are exponentiated.

There are a large number of possible terms that can be represented, however the most common are nearest-neighbour and one-local terms. Both cases are represented using the LocalOp structure:

julia> X,Y,Z = PAULI # A convenience `Tuple` containing the 3 Pauli matrices.(ComplexF64[0.0 + 0.0im 1.0 + 0.0im; 1.0 + 0.0im 0.0 + 0.0im], ComplexF64[0.0 + 0.0im 0.0 - 1.0im; 0.0 + 1.0im 0.0 + 0.0im], ComplexF64[1.0 + 0.0im 0.0 + 0.0im; 0.0 + 0.0im -1.0 + 0.0im])
julia> X0 = LocalOp(X) # One-local termLocalOp{1}:
 1.0 × | 2×2 Matrix{ComplexF64}: |
       |  0.0+0.0im  1.0+0.0im   |
       |  1.0+0.0im  0.0+0.0im   |
julia> ZZ = LocalOp(Z, Z) # Nearest neighbour termLocalOp{2}:
 1.0 × | 2×2 Matrix{ComplexF64}: | ⊗ | 2×2 Matrix{ComplexF64}: |
       |  1.0+0.0im   0.0+0.0im  |   |  1.0+0.0im   0.0+0.0im  |
       |  0.0+0.0im  -1.0+0.0im  |   |  0.0+0.0im  -1.0+0.0im  |

The terms defined above are objects in the following type hierarchy:

julia> typeof(X0) <: LocalOp{1} <: AbstractSingleTerm <: AbstractTermtrue
julia> typeof(ZZ) <: LocalOp{2} <: AbstractSingleTerm <: AbstractTermtrue

Terms can be scaled,

julia> 0.5 * ZZLocalOp{2}:
 0.5 × | 2×2 Matrix{ComplexF64}: | ⊗ | 2×2 Matrix{ComplexF64}: |
       |  1.0+0.0im   0.0+0.0im  |   |  1.0+0.0im   0.0+0.0im  |
       |  0.0+0.0im  -1.0+0.0im  |   |  0.0+0.0im  -1.0+0.0im  |

and added together to create Hamiltonian objects:

julia> H = 1.5 / 4 * ZZ - X0Hamiltonian containing:
  1 on-site term(s)
  1 local term(s)

which can also be scaled added to other terms

julia> H = 0.5 * H + LocalOp(X,Y) + LocalOp(Y,X)Hamiltonian containing:
  1 on-site term(s)
  3 local term(s)

Adding a Dissipator to a Hamiltonian constructs a Liouvillian object:

julia> L = H + Dissipator(Z, 0.5 * (X - im * Y))Liouvillian containing:
  1 on-site term(s)
  3 local term(s)
  2 lindblad operator(s)

A Dissipator can be constructed with any number of Lindblad jump operators as positional arguments to the constructor.

The Model container

Optionally, a Model object can be constructed which allows greater control over the Suzuki-Trotter decomposition of the time-evolution operator. Layers can be added to a Model both at time of construction and later using push!.

julia> XX = LocalOp(X,X)LocalOp{2}:
 1.0 × | 2×2 Matrix{ComplexF64}: | ⊗ | 2×2 Matrix{ComplexF64}: |
       |  0.0+0.0im  1.0+0.0im   |   |  0.0+0.0im  1.0+0.0im   |
       |  1.0+0.0im  0.0+0.0im   |   |  1.0+0.0im  0.0+0.0im   |
julia> YY = LocalOp(Y,Y)LocalOp{2}:
 1.0 × | 2×2 Matrix{ComplexF64}: | ⊗ | 2×2 Matrix{ComplexF64}: |
       |  0.0+0.0im  0.0-1.0im   |   |  0.0+0.0im  0.0-1.0im   |
       |  0.0+1.0im  0.0+0.0im   |   |  0.0+1.0im  0.0+0.0im   |
julia> model = Model(XX, YY)Model{Hamiltonian} with 2 Trotter layers
julia> push!(model, ZZ) # ZZ defined in the previous sectionModel{Hamiltonian} with 3 Trotter layers

If the second argument to push! is a layer of purely one-local terms (including terms of type Dissipator as TimeEvolutionPEPO.jl currently supports on-site dissipation only), then this layer is distributed evenly amongst all existing layers in model, multiplied by the appropriate prefactor. If this is undesired, then either add the one-local terms to the model first, or explicitly construct a Hamiltonian or Liouvillian containing these local terms before pushing this to model.

julia> model = Model(XX, YY)Model{Hamiltonian} with 2 Trotter layers
julia> push!(model, LocalOp(X))Model{Hamiltonian} with 2 Trotter layers
julia> push!(model, zero(Hamiltonian) + LocalOp(X))Model{Hamiltonian} with 3 Trotter layers
julia> push!(model, Hamiltonian(LocalOp(X)))Model{Hamiltonian} with 4 Trotter layers

Note, if one wishes to eventually push! a dissipator to an object of type Model, then it is best to explicitly specify the type parameter at construction as pushing a Dissipator or Liouvillian onto a Model{Hamiltonian} will result in an error:

julia> typeof(model)Model{Hamiltonian}
julia> push!(model, Dissipator(Z))ERROR: MethodError: no method matching Hamiltonian(::Liouvillian{Tuple{LocalOp{2}}}, ::Vector{LocalOp{1}})
The type `Hamiltonian` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  Hamiltonian(::Any)
   @ TimeEvolutionPEPO ~/work/TimeEvolutionPEPO.jl/TimeEvolutionPEPO.jl/src/model/hamiltonian.jl:14
  Hamiltonian(::T, ::Vector{LocalOp{1}}) where T<:Tuple{Vararg{AbstractSingleTerm}}
   @ TimeEvolutionPEPO ~/work/TimeEvolutionPEPO.jl/TimeEvolutionPEPO.jl/src/model/hamiltonian.jl:8
  Hamiltonian(::AbstractSingleTerm, ::Any)
   @ TimeEvolutionPEPO ~/work/TimeEvolutionPEPO.jl/TimeEvolutionPEPO.jl/src/model/hamiltonian.jl:12
  ...

As an example, consider the spin-$1/2$ Heisenberg Hamiltonian. Lets define a simple helper function to construct this model:

function heisenberg(; Jx, Jy, Jz)
    X,Y,Z = PAULI

    # Initialise an instance of the model object
    model = Model{Hamiltonian}()

    push!(model, Jx * LocalOp(X,X))
    push!(model, Jy * LocalOp(Y,Y))
    push!(model, Jz * LocalOp(Z,Z))

    return model
end

heisenberg (generic function with 1 method)

Now notice the difference between the model returned by out heisenberg function we just defined

julia> model_A = heisenberg(; Jx = 1, Jy = 1, Jz = 1)Model{Hamiltonian} with 3 Trotter layers

and the same model constructed by summing terms together with +:

julia> model_B = Model(XX + YY + ZZ)Model{Hamiltonian} with 1 Trotter layer

The former has 3 layers, whereas the latter only has 1! The difference manifests in the Trotter decomposition: model_A will result in the time evolution operator:

\[e^{\tau H_{A}} \approx e^{\tau X \otimes X} e^{\tau Y \otimes Y} e^{\tau Z \otimes Z}\]

whereas model_B:

\[e^{\tau H_{B}} \approx e^{\tau \left( X \otimes X + Y \otimes Y + Z \otimes Z \right)}\]

which can be important when individual Trotter layers have lower rank than the sum of the operators. Many methods can take advantage of this and construct operators of lower bond dimension, at the cost of more operators (and thus more truncation steps) in the Trotter decomposition. For example, for the TEBD method:

julia> state = PEPO(rand, 2, 1) # Random product state with D = 11×1 PEPO{Square, TensorKit.TensorMap{ComplexF64, TensorKit.ComplexSpace, 2, 4, Vector{ComplexF64}}, TensorKit.TensorMap{Float64, TensorKit.ComplexSpace, 1, 1, Vector{Float64}}}:
 TensorMap((ℂ^2 ⊗ (ℂ^2)') ← (ℂ^1 ⊗ ℂ^1 ⊗ (ℂ^1)' ⊗ (ℂ^1)')):
[:, :, 1, 1, 1, 1] =
 0.31897468754130676 + 0.7960962437497708im   …  0.48032278996161315 + 0.52036153041591im
   0.687919377972302 + 0.03204981420712483im      0.6286305893407358 + 0.7817236599513337im
julia> method = TEBD(compresstol = 10 * eps())Time-evolving block decimation (TEBD)
  compresstol: 2.220446049250313e-15
  alg: Simple update (SimpleUpdate)
    gauge: NoGauge()
    reduce: :pre
    zeta: Inf
  verbose: true
julia> Simulation(state, model_A, method; timestep = 0.1);[ Info: First-order trotter decomposition...
[ Info: Compressed to η=4 with error 1.1605376764834596e-15
[ Info: Compressed to η=4 with error 1.1839925039167004e-15
[ Info: Compressed to η=4 with error 7.27552540535717e-18
[ Info: Compressed to η=4 with error 1.1605376764834596e-15
[ Info: Compressed to η=4 with error 1.1839925039167004e-15
[ Info: Compressed to η=4 with error 7.27552540535717e-18
[ Info: ...of 6 operators
julia> Simulation(state, model_B, method; timestep = 0.1);[ Info: First-order trotter decomposition...
[ Info: Compressed to η=16 with error 0.0
[ Info: Compressed to η=16 with error 0.0
[ Info: ...of 2 operators

The option compresstol sets a threshold for truncating the bond dimension of an operator in the decomposition. One can see that the individually Trotterised operators in model_A essentially have a rank of 4, and therefore the bond dimension can be truncated down to 4 with no (numerical) error. In model_B, the sum of the terms has a much higher rank when exponentiated and cannot be truncated nearly as well.

Unfortunately, one-local terms and dissipation can lead to higher rank decompositions, even when operators are exponentiated separately.

julia> Simulation(state, Model(XX + Dissipator(Z)), method; timestep = 0.1);[ Info: First-order trotter decomposition...
[ Info: Compressed to η=8 with error 6.03108801515236e-17
[ Info: Compressed to η=8 with error 6.03108801515236e-17
[ Info: ...of 2 operators

To circumvent this, one can either put the one-local terms on a separate layer,

julia> Simulation(state, Model([XX, Dissipator(Z)]), method; timestep = 0.1);[ Info: First-order trotter decomposition...
[ Info: Compressed to η=4 with error 2.1152739990922743e-16
[ Info: Compressed to η=1 with error 8.94073142736696e-17
[ Info: Compressed to η=4 with error 2.1152739990922743e-16
[ Info: Compressed to η=1 with error 8.94073142736696e-17
[ Info: ...of 4 operators

or use one of the specific optimised terms as described in the next section.

Specific model components

In some cases, terms can have analytical expressions when exponentiated, benefit from more bespoke optimisation beyond simple low-rank compression, or are can be put in forms with specific properties that may be desirable. The model interface allows for this. For example, an Ising term coupled with optional local operators and a optional dissipation is specified with the Ising type, which defines a specific, optimised exponentiation which can be used instead of the more general LocalOp constructor. When usingIsing, low-rank forms of the Ising-type operators are obtained even when dissipation or one-local terms are present.

julia> model = Model(Ising(X, Dissipator(Z)))ERROR: MethodError: no method matching Ising(::Matrix{ComplexF64}, ::Vector{Any}, ::Dissipator)
The type `Ising` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  Ising(::LocalOp{2}, ::Any, ::Any)
   @ TimeEvolutionPEPO ~/work/TimeEvolutionPEPO.jl/TimeEvolutionPEPO.jl/src/model/ising.jl:40
  Ising(::Any, ::Any)
   @ TimeEvolutionPEPO ~/work/TimeEvolutionPEPO.jl/TimeEvolutionPEPO.jl/src/model/ising.jl:63
  Ising(::LocalOp{2}, ::Any, ::Nothing)
   @ TimeEvolutionPEPO ~/work/TimeEvolutionPEPO.jl/TimeEvolutionPEPO.jl/src/model/ising.jl:35
  ...
julia> Simulation(state, model, method; timestep = 0.1);[ Info: First-order trotter decomposition...
[ Info: Compressed to η=6 with error 2.4366760516432086e-15
[ Info: Compressed to η=16 with error 0.0
[ Info: Compressed to η=1 with error 2.350926934698294e-15
[ Info: Compressed to η=1 with error 2.350926934698294e-15
[ Info: Compressed to η=6 with error 2.4366760516432086e-15
[ Info: Compressed to η=16 with error 0.0
[ Info: Compressed to η=1 with error 2.350926934698294e-15
[ Info: Compressed to η=1 with error 2.350926934698294e-15
[ Info: ...of 8 operators