GKSL master equation

The GKSL master equation is general continuous-time open quantum evolution. The master equation base on to form of operators: the Hamiltonian, which describes the evolution of a closed system, and the Lindbladian operators which describes the evolution of an open system. Basic facts in the context of quantum stochastic walks can found in this preprint or its published version. For local interaction, where each Lindblad operator is a matrix with a single nonzero element, we provide a local_lind function which splits the matrix into mentioned operators.

Below we present documentation for essential functions used for simulating GKSL master equation.

QSWalk.bra
QSWalk.evolve
QSWalk.evolve_generator
QSWalk.evolve_operator
QSWalk.ket
QSWalk.ketbra
QSWalk.local_lind
QSWalk.proj
QSWalk.res
QSWalk.unres

Full docs

QSWalk.ket — Function

ket(index, size)

Return index-th base (column) vector in the size-dimensional vector space. To be consistent with Julia indexing, index = 1, 2, ..., size.

Examples

julia> ket(1, 2)
2-element SparseArrays.SparseVector{Int64,Int64} with 1 stored entry:
  [1]  =  1

source

QSWalk.bra — Function

bra(index, size)

Return index-th base row vector in the size-dimensional vector space, with index = 1, 2, ..., size.

Examples

julia> bra(1, 2)
1×2 LinearAlgebra.Adjoint{Int64,SparseArrays.SparseVector{Int64,Int64}}:
 1  0

source

QSWalk.ketbra — Function

ketbra(irow, icol, size)

Return matrix acting on size-dimensional vector space, irow, icol = 1, 2, ..., size. The matrix consists of single non-zero element equal to one, located at position (irow, icol).

Examples

julia> ketbra(1, 2, 3)
3×3 SparseArrays.SparseMatrixCSC{Int64,Int64} with 1 stored entry:
  [1, 2]  =  1

source

QSWalk.proj — Function

proj(index, size)

Return projector onto index-th base vector in size-dimensional vector space, with index = 1, 2, ..., size. This is equivalent to ketbra(index, index, size).

julia> proj(1, 2)
2×2 SparseArrays.SparseMatrixCSC{Int64,Int64} with 1 stored entry:
  [1, 1]  =  1

source

proj(vector)

Return projector onto the subspace spanned by vector vector.

Examples

julia> v = 1/sqrt(2) * (ket(1, 3)+ket(3, 3))
3-element SparseArrays.SparseVector{Float64,Int64} with 2 stored entries:
  [1]  =  0.707107
  [3]  =  0.707107

julia> proj(v)
3×3 SparseArrays.SparseMatrixCSC{Float64,Int64} with 4 stored entries:
  [1, 1]  =  0.5
  [3, 1]  =  0.5
  [1, 3]  =  0.5
  [3, 3]  =  0.5

source

QSWalk.res — Function

res(matrix)

Return vectorization of the matrix in the row order. This is equivalent to Base.vec(transpose(matrix)).

Examples

julia> M = reshape(1:9, (3, 3))'*1im
3×3 Array{Complex{Int64},2}:
 0+1im  0+2im  0+3im
 0+4im  0+5im  0+6im
 0+7im  0+8im  0+9im

julia> v = res(M)
9-element reshape(::LinearAlgebra.Transpose{Complex{Int64},Array{Complex{Int64},2}}, 9) with eltype Complex{Int64}:
 0 + 1im
 0 + 2im
 0 + 3im
 0 + 4im
 0 + 5im
 0 + 6im
 0 + 7im
 0 + 8im
 0 + 9im

julia> res(unres(v)) == v
true

source

QSWalk.unres — Function

unres(vector)

Return square matrix with elements from vector. The vector is expected to have perfect square number of arguments.

Examples

julia> v = collect(1:9)*im
9-element Array{Complex{Int64},1}:
 0 + 1im
 0 + 2im
 0 + 3im
 0 + 4im
 0 + 5im
 0 + 6im
 0 + 7im
 0 + 8im
 0 + 9im

julia> unres(v)
3×3 LinearAlgebra.Transpose{Complex{Int64},Array{Complex{Int64},2}}:
 0+1im  0+2im  0+3im
 0+4im  0+5im  0+6im
 0+7im  0+8im  0+9im

julia> res(unres(v)) == v
true

source

QSWalk.evolve_generator — Method

evolve_generator(H, L[, localH][, ω])

Create the generator for the evolution superoperator. Given Hamiltonian H, collection of Lindblad operator L, local Hamiltonian localH and scaling parameter ω, the generator is obtained as a sum

$-i(1-ω) (H ⊗ 1 - 1 ⊗ H) + ω (-i(localH ⊗ 1 - 1 ⊗ localH) + ∑(L ⊗ L̄ - 1/2(L^†L ⊗ 1 + 1 ⊗ L^T L̄ )))$

The last two arguments are optional.

If localH is not given, it defaults to sparse zero matrix of the size of H.

If ω is not given, both parts are taken with the same intensity and the global operator takes the form

$-i(H ⊗ 1 - 1 ⊗ H) + (-i(localH ⊗ 1 - 1 ⊗ localH) + ∑(L ⊗ L̄ - 1/2(L^†L ⊗ 1 + 1 ⊗ L^T L̄ )))$

Arguments

H: Hamiltonian, must be hermitian,
L: collection of Lindblad operators, each must be of the same size as H,
localH: local Hamiltonian, suggested for nonmoralized QS walk, must be hermitian and of the size of H,
ω: scaling parameter, should be in [0, 1].

Return

The generator matrix, which can be used in evolve function.

Examples

julia> H, L, localH = [0 1+im; 1-im 0], [0. 1; 0 0], [1. 0.; 0. 1.]
(Complex{Int64}[0 + 0im 1 + 1im; 1 - 1im 0 + 0im], [0.0 1.0; 0.0 0.0], [1.0 0.0; 0.0 1.0])

julia> evolve_generator(H, [L], localH, 1/2)
4×4 Array{Complex{Float64},2}:
  0.0+0.0im    0.5+0.5im    0.5-0.5im   0.5+0.0im
 -0.5+0.5im  -0.25+0.0im    0.0+0.0im   0.5-0.5im
 -0.5-0.5im    0.0+0.0im  -0.25+0.0im   0.5+0.5im
  0.0+0.0im   -0.5-0.5im   -0.5+0.5im  -0.5+0.0im

source

QSWalk.evolve_operator — Function

evolve_operator(evo_gen, time)

Return an exponent of time×evo_gen. This function is useful in the case of multiple initial states and fixed evo_gen and time, as it is faster to compute the exponent once and use it for evolving on different initial states.

Note: Parameter evo_gen must be of type Matrix. For type SparseMatrixCSC case different numerical approach is used. See function epmv in package Expokit.

Examples

julia> H, L = [0 1; 1 0], [[0 1; 0 0], [0 0; 1 0]]
([0 1; 1 0], Array{Int64,2}[[0 1; 0 0], [0 0; 1 0]])

julia> evolve_operator(evolve_generator(H, L), 4.0)
4×4 Array{Complex{Float64},2}:
 0.499815+0.0im                0.0+0.00127256im  …  0.500185+0.0im
      0.0+0.00127256im  0.00960957+0.0im                 0.0-0.00127256im
      0.0-0.00127256im  0.00870607+0.0im                 0.0+0.00127256im
 0.500185+0.0im                0.0-0.00127256im     0.499815+0.0im

source

QSWalk.evolve — Method

evolve(evogen, state, time)
evolve(evogen, state, tpoints)
evolve(evosuper, state)

Simulate the GKSL master equation according to the equation

$|result⟩⟩ = exp(time*evogen)|state⟩⟩$

where $|⋅⟩⟩$ denotes the vectorization.

Note: The function returns unvectorized result.

The evolution can be calculated using three different approaches.

In the simplest case the function accepts matrix evogen specifying the generator of the evolution, state describing the starting point of the evolution, and time specifying the time of the evolution.

Note: If evogen is of type Matrix, the exponent is calculated using exp function. If evogen is of type SparseMatrixCSC, expmv from Expokit.jl is used.

Alternatively, a list of point of time (tpoints) can be given. Points of time needs to be non-negative (you cannot go back in time). In this case a list of resulting states is returned.

Note: It is up to the user to provide evogen and state fulfilling the appropriate conditions. For the procedure to work correctly evogen should be generated by evolve_generator function and state should be a proper density matrix.

The third approach can be used if the superoperator is known. In this case argument evosuper can be specified. This argument can be generated by evolve_operator function. This is useful to simulate a fixed model of evolution in the case of multiple initial states and the same time point.

Examples

julia> H, L = [0 1; 1 0], [[0 1; 0 0], [0 0; 1 0]]
([0 1; 1 0], Array{Int64,2}[[0 1; 0 0], [0 0; 1 0]])

julia> Matrix(evolve(evolve_generator(H, L), proj(1, 2), 4.))
2×2 Array{Complex{Float64},2}:
 0.499815+0.0im              0.0+0.00127256im
      0.0-0.00127256im  0.500185+0.0im

julia> Matrix.(evolve(evolve_generator(H, L), proj(1, 2), [1., 2., 3., 4.]))
4-element Array{Array{Complex{Float64},2},1}:
 [0.43320327312718954 + 0.0im 0.0 + 0.10760477583156115im; 0.0 - 0.10760477583156117im 0.5667967268728107 + 0.0im]
 [0.48576602987832085 + 0.0im 0.0 - 0.017171799503962696im; 0.0 + 0.017171799503962682im 0.5142339701216796 + 0.0im]
 [0.5055971005015641 + 0.0im 0.0 - 0.002617013961759585im; 0.0 + 0.002617013961759615im 0.49440289949843613 + 0.0im]
 [0.49981547041444147 + 0.0im 0.0 + 0.0012725622225037549im; 0.0 - 0.0012725622225037724im 0.5001845295855593 + 0.0im]

julia> ev_op = evolve_operator(evolve_generator(H, L), 4.)
4×4 Array{Complex{Float64},2}:
 0.499815+0.0im                0.0+0.00127256im  …  0.500185+0.0im
      0.0+0.00127256im  0.00960957+0.0im                 0.0-0.00127256im
      0.0-0.00127256im  0.00870607+0.0im                 0.0+0.00127256im
 0.500185+0.0im                0.0-0.00127256im     0.499815+0.0im

julia> Matrix(evolve(ev_op, proj(1, 2)))
2×2 Array{Complex{Float64},2}:
 0.499815+0.0im              0.0+0.00127256im
      0.0-0.00127256im  0.500185+0.0im

source

QSWalk.local_lind — Function

local_lind(A[; epsilon])

Split the elements of the matrix A into a collection of sparse matrices with exactly one non-zero element. Martices are created if the absolute value of the nonzero element is there are not smaller than epsilon, where epsilon should be nonnegative. The epsilon defaults to eps() if not specified.

Examples

julia> A = [1. 2.; 3. 4.]
2×2 Array{Float64,2}:
 1.0  2.0
 3.0  4.0

julia> local_lind(A)
4-element Array{SparseArrays.SparseMatrixCSC{Float64,Ti} where Ti<:Integer,1}:

  [1, 1]  =  1.0

  [1, 2]  =  2.0

  [2, 1]  =  3.0

  [2, 2]  =  4.0

julia> local_lind(A, epsilon = 1.5)
3-element Array{SparseArrays.SparseMatrixCSC{Float64,Ti} where Ti<:Integer,1}:

  [1, 2]  =  2.0

  [2, 1]  =  3.0

  [2, 2]  =  4.0

source