API

BraQuantumObject <: QuantumObjectType

Constructor representing a bra state $\langle\psi|$.

const Bra = BraQuantumObject()

A constant representing the type of BraQuantumObject: a bra state $\langle\psi|$

KetQuantumObject <: QuantumObjectType

Constructor representing a ket state $|\psi\rangle$.

const Ket = KetQuantumObject()

A constant representing the type of KetQuantumObject: a ket state $|\psi\rangle$

OperatorQuantumObject <: QuantumObjectType

Constructor representing an operator $\hat{O}$.

const Operator = OperatorQuantumObject()

A constant representing the type of OperatorQuantumObject: an operator $\hat{O}$

OperatorBraQuantumObject <: QuantumObjectType

Constructor representing a bra state in the SuperOperator formalism $\langle\langle\rho|$.

const OperatorBra = OperatorBraQuantumObject()

A constant representing the type of OperatorBraQuantumObject: a bra state in the SuperOperator formalism $\langle\langle\rho|$.

OperatorKetQuantumObject <: QuantumObjectType

Constructor representing a ket state in the SuperOperator formalism $|\rho\rangle\rangle$.

const OperatorKet = OperatorKetQuantumObject()

A constant representing the type of OperatorKetQuantumObject: a ket state in the SuperOperator formalism $|\rho\rangle\rangle$

SuperOperatorQuantumObject <: QuantumObjectType

Constructor representing a super-operator $\hat{\mathcal{O}}$ acting on vectorized density operator matrices.

const SuperOperator = SuperOperatorQuantumObject()

A constant representing the type of SuperOperatorQuantumObject: a super-operator $\hat{\mathcal{O}}$ acting on vectorized density operator matrices

struct QuantumObject{MT<:AbstractArray,ObjType<:QuantumObjectType,N}
    data::MT
    type::ObjType
    dims::SVector{N, Int}
end

Julia struct representing any quantum objects.

Examples

julia> a = destroy(20)
Quantum Object:   type=Operator   dims=[20]   size=(20, 20)   ishermitian=false
20×20 SparseMatrixCSC{ComplexF64, Int64} with 19 stored entries:
⎡⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⎦

julia> a isa QuantumObject
true

struct OperatorSum

A constructor to represent a sum of operators $\sum_i c_i \hat{O}_i$ with a list of coefficients $c_i$ and a list of operators $\hat{O}_i$.

This is very useful when we have to update only the coefficients, without allocating memory by performing the sum of the operators.

size(A::QuantumObject)
size(A::QuantumObject, idx::Int)

Returns a tuple containing each dimensions of the array in the QuantumObject.

Optionally, you can specify an index (idx) to just get the corresponding dimension of the array.

eltype(A::QuantumObject)

Returns the elements type of the matrix or vector corresponding to the QuantumObject A.

length(A::QuantumObject)

Returns the length of the matrix or vector corresponding to the QuantumObject A.

isbra(A::QuantumObject)

Checks if the QuantumObject A is a BraQuantumObject.

isket(A::QuantumObject)

Checks if the QuantumObject A is a KetQuantumObject.

isoper(A::QuantumObject)

Checks if the QuantumObject A is a OperatorQuantumObject.

isoperbra(A::QuantumObject)

Checks if the QuantumObject A is a OperatorBraQuantumObject.

isoperket(A::QuantumObject)

Checks if the QuantumObject A is a OperatorKetQuantumObject.

issuper(A::QuantumObject)

Checks if the QuantumObject A is a SuperOperatorQuantumObject.

ishermitian(A::QuantumObject)

Test whether the QuantumObject is Hermitian.

issymmetric(A::QuantumObject)

Test whether the QuantumObject is symmetric.

isposdef(A::QuantumObject)

Test whether the QuantumObject is positive definite (and Hermitian) by trying to perform a Cholesky factorization of A.

isunitary(U::QuantumObject; kwargs...)

Test whether the QuantumObject $U$ is unitary operator. This function calls Base.isapprox to test whether $U U^\dagger$ is approximately equal to identity operator.

Note that all the keyword arguments will be passed to Base.isapprox.

conj(A::QuantumObject)

Return the element-wise complex conjugation of the QuantumObject.

transpose(A::QuantumObject)

Lazy matrix transpose of the QuantumObject.

A'
adjoint(A::QuantumObject)

Lazy adjoint (conjugate transposition) of the QuantumObject

Note that A' is a synonym for adjoint(A)

dot(A::QuantumObject, B::QuantumObject)

Compute the dot product between two QuantumObject: $\langle A | B \rangle$

Note that A and B should be Ket or OperatorKet

A ⋅ B (where ⋅ can be typed by tab-completing \cdot in the REPL) is a synonym for dot(A, B)

dot(i::QuantumObject, A::QuantumObject j::QuantumObject)

Compute the generalized dot product dot(i, A*j) between three QuantumObject: $\langle i | \hat{A} | j \rangle$

Supports the following inputs:

• A is in the type of Operator, with i and j are both Ket.
• A is in the type of SuperOperator, with i and j are both OperatorKet

√(A)
sqrt(A::QuantumObject)

Matrix square root of QuantumObject

Note that √(A) is a synonym for sqrt(A)

log(A::QuantumObject)

Matrix logarithm of QuantumObject

Note that this function only supports for Operator and SuperOperator

exp(A::QuantumObject)

Matrix exponential of QuantumObject

Note that this function only supports for Operator and SuperOperator

sin(A::QuantumObject)

Matrix sine of QuantumObject, defined as

$\sin \left( \hat{A} \right) = \frac{e^{i \hat{A}} - e^{-i \hat{A}}}{2 i}$

Note that this function only supports for Operator and SuperOperator

cos(A::QuantumObject)

Matrix cosine of QuantumObject, defined as

$\cos \left( \hat{A} \right) = \frac{e^{i \hat{A}} + e^{-i \hat{A}}}{2}$

Note that this function only supports for Operator and SuperOperator

tr(A::QuantumObject)

Returns the trace of QuantumObject.

Note that this function only supports for Operator and SuperOperator

Examples

julia> a = destroy(20)
Quantum Object:   type=Operator   dims=[20]   size=(20, 20)   ishermitian=false
20×20 SparseMatrixCSC{ComplexF64, Int64} with 19 stored entries:
⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢

julia> tr(a' * a)
190.0 + 0.0im

svdvals(A::QuantumObject)

Return the singular values of a QuantumObject in descending order

norm(A::QuantumObject, p::Real)

Return the standard vector p-norm or Schatten p-norm of a QuantumObject depending on the type of A:

• If A is either Ket, Bra, OperatorKet, or OperatorBra, returns the standard vector p-norm (default p=2) of A.
• If A is either Operator or SuperOperator, returns Schatten p-norm (default p=1) of A.

Examples

julia> ψ = fock(10, 2)
Quantum Object:   type=Ket   dims=[10]   size=(10,)
10-element Vector{ComplexF64}:
 0.0 + 0.0im
 0.0 + 0.0im
 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

julia> norm(ψ)
1.0

normalize(A::QuantumObject, p::Real)

Return normalized QuantumObject so that its p-norm equals to unity, i.e. norm(A, p) == 1.

Support for the following types of QuantumObject:

• If A is Ket or Bra, default p = 2
• If A is Operator, default p = 1

Also, see norm about its definition for different types of QuantumObject.

normalize!(A::QuantumObject, p::Real)

Normalize QuantumObject in-place so that its p-norm equals to unity, i.e. norm(A, p) == 1.

Support for the following types of QuantumObject:

• If A is Ket or Bra, default p = 2
• If A is Operator, default p = 1

Also, see norm about its definition for different types of QuantumObject.

inv(A::QuantumObject)

Matrix inverse of the QuantumObject

diag(A::QuantumObject, k::Int=0)

Return the k-th diagonal elements of a matrix-type QuantumObject

Note that this function only supports for Operator and SuperOperator

proj(ψ::QuantumObject)

Return the projector for a Ket or Bra type of QuantumObject

ptrace(QO::QuantumObject, sel)

Partial trace of a quantum state QO leaving only the dimensions with the indices present in the sel vector.

Note that this function will always return Operator. No matter the input QuantumObject is a Ket, Bra, or Operator.

Examples

Two qubits in the state $\ket{\psi} = \ket{e,g}$:

julia> ψ = kron(fock(2,0), fock(2,1))
Quantum Object:   type=Ket   dims=[2, 2]   size=(4,)
4-element Vector{ComplexF64}:
 0.0 + 0.0im
 1.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im

julia> ptrace(ψ, 2)
Quantum Object:   type=Operator   dims=[2]   size=(2, 2)   ishermitian=true
2×2 Matrix{ComplexF64}:
 0.0+0.0im  0.0+0.0im
 0.0+0.0im  1.0+0.0im

or in an entangled state $\ket{\psi} = \frac{1}{\sqrt{2}} \left( \ket{e,e} + \ket{g,g} \right)$:

julia> ψ = 1 / √2 * (kron(fock(2,0), fock(2,0)) + kron(fock(2,1), fock(2,1)))
Quantum Object:   type=Ket   dims=[2, 2]   size=(4,)
4-element Vector{ComplexF64}:
 0.7071067811865475 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
 0.7071067811865475 + 0.0im

julia> ptrace(ψ, 1)
Quantum Object:   type=Operator   dims=[2]   size=(2, 2)   ishermitian=true
2×2 Matrix{ComplexF64}:
 0.5+0.0im  0.0+0.0im
 0.0+0.0im  0.5+0.0im

purity(ρ::QuantumObject)

Calculate the purity of a QuantumObject: $\textrm{Tr}(\rho^2)$

Note that this function only supports for Ket, Bra, and Operator

permute(A::QuantumObject, order::Union{AbstractVector{Int},Tuple})

Permute the tensor structure of a QuantumObject A according to the specified order list

Note that this method currently works for Ket, Bra, and Operator types of QuantumObject.

Examples

If order = [2, 1, 3], the Hilbert space structure will be re-arranged: $\mathcal{H}_1 \otimes \mathcal{H}_2 \otimes \mathcal{H}_3 \rightarrow \mathcal{H}_2 \otimes \mathcal{H}_1 \otimes \mathcal{H}_3$.

julia> ψ1 = fock(2, 0)
julia> ψ2 = fock(3, 1)
julia> ψ3 = fock(4, 2)
julia> ψ_123 = tensor(ψ1, ψ2, ψ3)
julia> permute(ψ_123, [2, 1, 3]) ≈ tensor(ψ2, ψ1, ψ3)
true

tidyup(A::QuantumObject, tol::Real=1e-14)

Given a QuantumObject A, check the real and imaginary parts of each element separately. Remove the real or imaginary value if its absolute value is less than tol.

tidyup!(A::QuantumObject, tol::Real=1e-14)

Given a QuantumObject A, check the real and imaginary parts of each element separately. Remove the real or imaginary value if its absolute value is less than tol.

Note that this function is an in-place version of tidyup.

get_data(A::QuantumObject)

Returns the data of a QuantumObject.

get_coherence(ψ::QuantumObject)

Get the coherence value $\alpha$ by measuring the expectation value of the destruction operator $\hat{a}$ on a state $\ket{\psi}$ or a density matrix $\hat{\rho}$.

It returns both $\alpha$ and the corresponding state with the coherence removed: $\ket{\delta_\alpha} = \exp ( \alpha^* \hat{a} - \alpha \hat{a}^\dagger ) \ket{\psi}$ for a pure state, and $\hat{\rho_\alpha} = \exp ( \alpha^* \hat{a} - \alpha \hat{a}^\dagger ) \hat{\rho} \exp ( -\bar{\alpha} \hat{a} + \alpha \hat{a}^\dagger )$ for a density matrix. These states correspond to the quantum fluctuations around the coherent state $\ket{\alpha}$ or $|\alpha\rangle\langle\alpha|$.

partial_transpose(ρ::QuantumObject, mask::Vector{Bool})

Return the partial transpose of a density matrix $\rho$, where mask is an array/vector with length that equals the length of ρ.dims. The elements in mask are boolean (true or false) which indicates whether or not the corresponding subsystem should be transposed.

Arguments

• ρ::QuantumObject: The density matrix (ρ.type must be OperatorQuantumObject).
• mask::Vector{Bool}: A boolean vector selects which subsystems should be transposed.

Returns

• ρ_pt::QuantumObject: The density matrix with the selected subsystems transposed.

struct EigsolveResult{T1<:Vector{<:Number}, T2<:AbstractMatrix{<:Number}, ObjType<:Union{Nothing,OperatorQuantumObject,SuperOperatorQuantumObject},N}
    values::T1
    vectors::T2
    type::ObjType
    dims::SVector{N,Int}
    iter::Int
    numops::Int
    converged::Bool
end

A struct containing the eigenvalues, the eigenvectors, and some information from the solver

Fields

• values::AbstractVector: the eigenvalues
• vectors::AbstractMatrix: the transformation matrix (eigenvectors)
• type::Union{Nothing,QuantumObjectType}: the type of QuantumObject, or nothing means solving eigen equation for general matrix
• dims::SVector: the dims of QuantumObject
• iter::Int: the number of iteration during the solving process
• numops::Int : number of times the linear map was applied in krylov methods
• converged::Bool: Whether the result is converged

Examples

One can obtain the eigenvalues and the corresponding QuantumObject-type eigenvectors by:

julia> result = eigenstates(sigmax());

julia> λ, ψ, T = result;

julia> λ
2-element Vector{ComplexF64}:
 -1.0 + 0.0im
  1.0 + 0.0im

julia> ψ
2-element Vector{QuantumObject{Vector{ComplexF64}, KetQuantumObject}}:
 QuantumObject{Vector{ComplexF64}, KetQuantumObject}(ComplexF64[-0.7071067811865475 + 0.0im, 0.7071067811865475 + 0.0im], KetQuantumObject(), [2])
 QuantumObject{Vector{ComplexF64}, KetQuantumObject}(ComplexF64[0.7071067811865475 + 0.0im, 0.7071067811865475 + 0.0im], KetQuantumObject(), [2])

julia> T
2×2 Matrix{ComplexF64}:
 -0.707107+0.0im  0.707107+0.0im
  0.707107+0.0im  0.707107+0.0im

eigenenergies(A::QuantumObject; sparse::Bool=false, kwargs...)

Calculate the eigenenergies

Arguments

• A::QuantumObject: the QuantumObject to solve eigenvalues
• sparse::Bool: if false call eigvals(A::QuantumObject; kwargs...), otherwise call eigsolve. Default to false.
• kwargs: Additional keyword arguments passed to the solver

Returns

• ::Vector{<:Number}: a list of eigenvalues

eigenstates(A::QuantumObject; sparse::Bool=false, kwargs...)

Calculate the eigenvalues and corresponding eigenvectors

Arguments

• A::QuantumObject: the QuantumObject to solve eigenvalues and eigenvectors
• sparse::Bool: if false call eigen(A::QuantumObject; kwargs...), otherwise call eigsolve. Default to false.
• kwargs: Additional keyword arguments passed to the solver

Returns

• ::EigsolveResult: containing the eigenvalues, the eigenvectors, and some information from the solver. see also EigsolveResult

LinearAlgebra.eigen(A::QuantumObject; kwargs...)

Calculates the eigenvalues and eigenvectors of the QuantumObject A using the Julia LinearAlgebra package.

julia> a = destroy(5);

julia> H = a + a'
Quantum Object:   type=Operator   dims=[5]   size=(5, 5)   ishermitian=true
5×5 SparseMatrixCSC{ComplexF64, Int64} with 8 stored entries:
     ⋅          1.0+0.0im          ⋅              ⋅          ⋅
 1.0+0.0im          ⋅      1.41421+0.0im          ⋅          ⋅
     ⋅      1.41421+0.0im          ⋅      1.73205+0.0im      ⋅
     ⋅              ⋅      1.73205+0.0im          ⋅      2.0+0.0im
     ⋅              ⋅              ⋅          2.0+0.0im      ⋅

julia> E, ψ, U = eigen(H)
EigsolveResult:   type=Operator   dims=[5]
values:
5-element Vector{Float64}:
 -2.8569700138728
 -1.3556261799742608
  1.3322676295501878e-15
  1.3556261799742677
  2.8569700138728056
vectors:
5×5 Matrix{ComplexF64}:
  0.106101+0.0im  -0.471249-0.0im  …   0.471249-0.0im  0.106101-0.0im
 -0.303127-0.0im   0.638838+0.0im      0.638838+0.0im  0.303127-0.0im
  0.537348+0.0im  -0.279149-0.0im      0.279149-0.0im  0.537348-0.0im
 -0.638838-0.0im  -0.303127-0.0im     -0.303127-0.0im  0.638838+0.0im
  0.447214+0.0im   0.447214+0.0im     -0.447214-0.0im  0.447214-0.0im

julia> expect(H, ψ[1]) ≈ E[1]
true

LinearAlgebra.eigvals(A::QuantumObject; kwargs...)

Same as eigen(A::QuantumObject; kwargs...) but for only the eigenvalues.

eigsolve(A::QuantumObject; 
    v0::Union{Nothing,AbstractVector}=nothing, 
    sigma::Union{Nothing, Real}=nothing,
    k::Int = 1,
    krylovdim::Int = max(20, 2*k+1),
    tol::Real = 1e-8,
    maxiter::Int = 200,
    solver::Union{Nothing, SciMLLinearSolveAlgorithm} = nothing,
    kwargs...)

Solve for the eigenvalues and eigenvectors of a matrix A using the Arnoldi method.

Notes

• For more details about solver and extra kwargs, please refer to LinearSolve.jl

Returns

• EigsolveResult: A struct containing the eigenvalues, the eigenvectors, and some information about the eigsolver

eigsolve_al(H::QuantumObject,
    T::Real, c_ops::Union{Nothing,AbstractVector}=nothing;
    alg::OrdinaryDiffEqAlgorithm=Tsit5(),
    H_t::Union{Nothing,Function}=nothing,
    params::NamedTuple=NamedTuple(),
    ρ0::Union{Nothing, AbstractMatrix} = nothing,
    k::Int=1,
    krylovdim::Int=min(10, size(H, 1)),
    maxiter::Int=200,
    eigstol::Real=1e-6,
    kwargs...)

Solve the eigenvalue problem for a Liouvillian superoperator L using the Arnoldi-Lindblad method.

Arguments

• H: The Hamiltonian (or directly the Liouvillian) of the system.
• T: The time at which to evaluate the time evolution
• c_ops: A vector of collapse operators. Default is nothing meaning the system is closed.
• alg: The differential equation solver algorithm
• H_t: A function H_t(t) that returns the additional term at time t
• params: A dictionary of additional parameters
• ρ0: The initial density matrix. If not specified, a random density matrix is used
• k: The number of eigenvalues to compute
• krylovdim: The dimension of the Krylov subspace
• maxiter: The maximum number of iterations for the eigsolver
• eigstol: The tolerance for the eigsolver
• kwargs: Additional keyword arguments passed to the differential equation solver

Notes

• For more details about alg please refer to DifferentialEquations.jl (ODE Solvers)
• For more details about kwargs please refer to DifferentialEquations.jl (Keyword Arguments)

Returns

• EigsolveResult: A struct containing the eigenvalues, the eigenvectors, and some information about the eigsolver

References

• [1] Minganti, F., & Huybrechts, D. (2022). Arnoldi-Lindblad time evolution: Faster-than-the-clock algorithm for the spectrum of time-independent and Floquet open quantum systems. Quantum, 6, 649.

ket2dm(ψ::QuantumObject)

Transform the ket state $\ket{\psi}$ into a pure density matrix $\hat{\rho} = \dyad{\psi}$.

expect(O::QuantumObject, ψ::Union{QuantumObject,Vector{QuantumObject}})

Expectation value of the Operator O with the state ψ. The state can be a Ket, Bra or Operator.

If ψ is a Ket or Bra, the function calculates $\langle\psi|\hat{O}|\psi\rangle$.

If ψ is a density matrix (Operator), the function calculates $\textrm{Tr} \left[ \hat{O} \hat{\psi} \right]$

The function returns a real number if O is of Hermitian type or Symmetric type, and returns a complex number otherwise. You can make an operator O hermitian by using Hermitian(O).

Note that ψ can also be given as a list of QuantumObject, it returns a list of expectation values.

Examples

julia> ψ = 1 / √2 * (fock(10,2) + fock(10,4));

julia> a = destroy(10);

julia> expect(a' * a, ψ) |> round
3.0 + 0.0im

julia> expect(Hermitian(a' * a), ψ) |> round
3.0

variance(O::QuantumObject, ψ::Union{QuantumObject,Vector{QuantumObject}})

Variance of the Operator O: $\langle\hat{O}^2\rangle - \langle\hat{O}\rangle^2$,

where $\langle\hat{O}\rangle$ is the expectation value of O with the state ψ (see also expect), and the state ψ can be a Ket, Bra or Operator.

The function returns a real number if O is hermitian, and returns a complex number otherwise.

Note that ψ can also be given as a list of QuantumObject, it returns a list of expectation values.

kron(A::QuantumObject, B::QuantumObject, ...)

Returns the Kronecker product $\hat{A} \otimes \hat{B} \otimes \cdots$.

Examples

julia> a = destroy(20)
Quantum Object:   type=Operator   dims=[20]   size=(20, 20)   ishermitian=false
20×20 SparseMatrixCSC{ComplexF64, Int64} with 19 stored entries:
⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢

julia> kron(a, a)
Quantum Object:   type=Operator   dims=[20, 20]   size=(400, 400)   ishermitian=false
400×400 SparseMatrixCSC{ComplexF64, Int64} with 361 stored entries:
⠀⠀⠘⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⠦

sparse_to_dense(A::QuantumObject)

Converts a sparse QuantumObject to a dense QuantumObject.

dense_to_sparse(A::QuantumObject)

Converts a dense QuantumObject to a sparse QuantumObject.

vec2mat(A::AbstractVector)

Converts a vector to a matrix.

vec2mat(A::QuantumObject)

Convert a quantum object from vector (OperatorKetQuantumObject-type) to matrix (OperatorQuantumObject-type)

mat2vec(A::QuantumObject)

Convert a quantum object from matrix (OperatorQuantumObject-type) to vector (OperatorKetQuantumObject-type)

mat2vec(A::AbstractMatrix)

Converts a matrix to a vector.

zero_ket(dimensions)

Returns a zero Ket vector with given argument dimensions.

The dimensions can be either the following types:

• dimensions::Int: Number of basis states in the Hilbert space.
• dimensions::Union{AbstractVector{Int}, Tuple}: list of dimensions representing the each number of basis in the subsystems.

fock(N::Int, j::Int=0; dims::Union{Int,AbstractVector{Int},Tuple}=N, sparse::Union{Bool,Val}=Val(false))

Generates a fock state $\ket{\psi}$ of dimension N.

It is also possible to specify the list of dimensions dims if different subsystems are present.

basis(N::Int, j::Int = 0; dims::Union{Int,AbstractVector{Int},Tuple}=N)

Generates a fock state like fock.

It is also possible to specify the list of dimensions dims if different subsystems are present.

coherent(N::Int, α::Number)

Generates a coherent state $|\alpha\rangle$, which is defined as an eigenvector of the bosonic annihilation operator $\hat{a} |\alpha\rangle = \alpha |\alpha\rangle$.

This state is constructed via the displacement operator displace and zero-fock state fock: $|\alpha\rangle = \hat{D}(\alpha) |0\rangle$

rand_ket(dimensions)

Generate a random normalized Ket vector with given argument dimensions.

The dimensions can be either the following types:

• dimensions::Int: Number of basis states in the Hilbert space.
• dimensions::Union{AbstractVector{Int},Tuple}: list of dimensions representing the each number of basis in the subsystems.

fock_dm(N::Int, j::Int=0; dims::Union{Int,AbstractVector{Int},Tuple}=N, sparse::Union{Bool,Val}=Val(false))

Density matrix representation of a Fock state.

Constructed via outer product of fock.

coherent_dm(N::Int, α::Number)

Density matrix representation of a coherent state.

Constructed via outer product of coherent.

thermal_dm(N::Int, n::Real; sparse::Union{Bool,Val}=Val(false))

Density matrix for a thermal state (generating thermal state probabilities) with the following arguments:

• N::Int: Number of basis states in the Hilbert space
• n::Real: Expectation value for number of particles in the thermal state.
• sparse::Union{Bool,Val}: If true, return a sparse matrix representation.

maximally_mixed_dm(dimensions)

Returns the maximally mixed density matrix with given argument dimensions.

The dimensions can be either the following types:

• dimensions::Int: Number of basis states in the Hilbert space.
• dimensions::Union{AbstractVector{Int},Tuple}: list of dimensions representing the each number of basis in the subsystems.

rand_dm(dimensions; rank::Int=prod(dimensions))

Generate a random density matrix from Ginibre ensemble with given argument dimensions and rank, ensuring that it is positive semi-definite and trace equals to 1.

The dimensions can be either the following types:

• dimensions::Int: Number of basis states in the Hilbert space.
• dimensions::Union{AbstractVector{Int},Tuple}: list of dimensions representing the each number of basis in the subsystems.

The default keyword argument rank = prod(dimensions) (full rank).

References

• J. Ginibre, Statistical ensembles of complex, quaternion, and real matrices, Journal of Mathematical Physics 6.3 (1965): 440-449
• K. Życzkowski, et al., Generating random density matrices, Journal of Mathematical Physics 52, 062201 (2011)

spin_state(j::Real, m::Real)

Generate the spin state: $|j, m\rangle$

The eigenstate of the Spin-j $\hat{S}_z$ operator with eigenvalue m, where where j is the spin quantum number and can be a non-negative integer or half-integer

See also jmat.

spin_coherent(j::Real, θ::Real, ϕ::Real)

Generate the coherent spin state (rotation of the $|j, j\rangle$ state), namely

\[|\theta, \phi \rangle = \hat{R}(\theta, \phi) |j, j\rangle\]

where the rotation operator is defined as

\[\hat{R}(\theta, \phi) = \exp \left( \frac{\theta}{2} (\hat{S}_- e^{i\phi} - \hat{S}_+ e^{-i\phi}) \right)\]

and $\hat{S}_\pm$ are plus and minus Spin-j operators, respectively.

Arguments

• j::Real: The spin quantum number and can be a non-negative integer or half-integer
• θ::Real: rotation angle from z-axis
• ϕ::Real: rotation angle from x-axis

See also jmat and spin_state.

Reference

• Robert Jones, Spin Coherent States and Statistical Physics

bell_state(x::Union{Int}, z::Union{Int})

Return the Bell state depending on the arguments (x, z):

• (0, 0): $| \Phi^+ \rangle = ( |00\rangle + |11\rangle ) / \sqrt{2}$
• (0, 1): $| \Phi^- \rangle = ( |00\rangle - |11\rangle ) / \sqrt{2}$
• (1, 0): $| \Psi^+ \rangle = ( |01\rangle + |10\rangle ) / \sqrt{2}$
• (1, 1): $| \Psi^- \rangle = ( |01\rangle - |10\rangle ) / \sqrt{2}$

Here, x = 1 (z = 1) means applying Pauli-$X$ ( Pauli-$Z$) unitary transformation on $| \Phi^+ \rangle$.

Example

julia> bell_state(0, 0)
Quantum Object:   type=Ket   dims=[2, 2]   size=(4,)
4-element Vector{ComplexF64}:
 0.7071067811865475 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
 0.7071067811865475 + 0.0im

julia> bell_state(Val(1), Val(0))
Quantum Object:   type=Ket   dims=[2, 2]   size=(4,)
4-element Vector{ComplexF64}:
                0.0 + 0.0im
 0.7071067811865475 + 0.0im
 0.7071067811865475 + 0.0im
                0.0 + 0.0im

singlet_state()

Return the two particle singlet state: $\frac{1}{\sqrt{2}} ( |01\rangle - |10\rangle )$

triplet_states()

Return a list of the two particle triplet states:

• $|11\rangle$
• $( |01\rangle + |10\rangle ) / \sqrt{2}$
• $|00\rangle$

w_state(n::Union{Int,Val})

Returns the n-qubit W-state:

\[\frac{1}{\sqrt{n}} \left( |100...0\rangle + |010...0\rangle + \cdots + |00...01\rangle \right)\]

ghz_state(n::Union{Int,Val}; d::Int=2)

Returns the generalized n-qudit Greenberger–Horne–Zeilinger (GHZ) state:

\[\frac{1}{\sqrt{d}} \sum_{i=0}^{d-1} | i \rangle \otimes \cdots \otimes | i \rangle\]

Here, d specifies the dimension of each qudit. Default to d=2 (qubit).

rand_unitary(dimensions, distribution=Val(:haar))

Returns a random unitary QuantumObject.

The dimensions can be either the following types:

• dimensions::Int: Number of basis states in the Hilbert space.
• dimensions::Union{AbstractVector{Int},Tuple}: list of dimensions representing the each number of basis in the subsystems.

The distribution specifies which of the method used to obtain the unitary matrix:

• :haar: Haar random unitary matrix using the algorithm from reference 1
• :exp: Uses $\exp(-i\hat{H})$, where $\hat{H}$ is a randomly generated Hermitian operator.

References

1. F. Mezzadri, How to generate random matrices from the classical compact groups, arXiv:math-ph/0609050 (2007)

sigmap()

Pauli ladder operator $\hat{\sigma}_+ = (\hat{\sigma}_x + i \hat{\sigma}_y) / 2$.

See also jmat.

sigmam()

Pauli ladder operator $\hat{\sigma}_- = (\hat{\sigma}_x - i \hat{\sigma}_y) / 2$.

See also jmat.

sigmax()

Pauli operator $\hat{\sigma}_x = \hat{\sigma}_- + \hat{\sigma}_+$.

See also jmat.

sigmay()

Pauli operator $\hat{\sigma}_y = i \left( \hat{\sigma}_- - \hat{\sigma}_+ \right)$.

See also jmat.

sigmaz()

Pauli operator $\hat{\sigma}_z = \comm{\hat{\sigma}_+}{\hat{\sigma}_-}$.

See also jmat.

jmat(j::Real, which::Union{Symbol,Val})

Generate higher-order Spin-j operators, where j is the spin quantum number and can be a non-negative integer or half-integer

The parameter which specifies which of the following operator to return.

• :x: $\hat{S}_x$
• :y: $\hat{S}_y$
• :z: $\hat{S}_z$
• :+: $\hat{S}_+$
• :-: $\hat{S}_-$

Note that if the parameter which is not specified, returns a set of Spin-j operators: $(\hat{S}_x, \hat{S}_y, \hat{S}_z)$

Examples

julia> jmat(0.5, :x)
Quantum Object:   type=Operator   dims=[2]   size=(2, 2)   ishermitian=true
2×2 SparseMatrixCSC{ComplexF64, Int64} with 2 stored entries:
     ⋅      0.5+0.0im
 0.5+0.0im      ⋅

julia> jmat(0.5, :-)
Quantum Object:   type=Operator   dims=[2]   size=(2, 2)   ishermitian=false
2×2 SparseMatrixCSC{ComplexF64, Int64} with 1 stored entry:
     ⋅          ⋅    
 1.0+0.0im      ⋅

julia> jmat(1.5, Val(:z))
Quantum Object:   type=Operator   dims=[4]   size=(4, 4)   ishermitian=true
4×4 SparseMatrixCSC{ComplexF64, Int64} with 4 stored entries:
 1.5+0.0im      ⋅           ⋅           ⋅    
     ⋅      0.5+0.0im       ⋅           ⋅    
     ⋅          ⋅      -0.5+0.0im       ⋅    
     ⋅          ⋅           ⋅      -1.5+0.0im

spin_Jx(j::Real)

$\hat{S}_x$ operator for Spin-j, where j is the spin quantum number and can be a non-negative integer or half-integer.

See also jmat.

spin_Jy(j::Real)

$\hat{S}_y$ operator for Spin-j, where j is the spin quantum number and can be a non-negative integer or half-integer.

See also jmat.

spin_Jz(j::Real)

$\hat{S}_z$ operator for Spin-j, where j is the spin quantum number and can be a non-negative integer or half-integer.

See also jmat.

spin_Jm(j::Real)

$\hat{S}_-$ operator for Spin-j, where j is the spin quantum number and can be a non-negative integer or half-integer.

See also jmat.

spin_Jp(j::Real)

$\hat{S}_+$ operator for Spin-j, where j is the spin quantum number and can be a non-negative integer or half-integer.

See also jmat.

spin_J_set(j::Real)

A set of Spin-j operators $(\hat{S}_x, \hat{S}_y, \hat{S}_z)$, where j is the spin quantum number and can be a non-negative integer or half-integer.

Note that this functions is same as jmat(j). See also jmat.

destroy(N::Int)

Bosonic annihilation operator with Hilbert space cutoff N.

This operator acts on a fock state as $\hat{a} \ket{n} = \sqrt{n} \ket{n-1}$.

Examples

julia> a = destroy(20)
Quantum Object:   type=Operator   dims=[20]   size=(20, 20)   ishermitian=false
20×20 SparseMatrixCSC{ComplexF64, Int64} with 19 stored entries:
⎡⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢⎦

julia> fock(20, 3)' * a * fock(20, 4)
2.0 + 0.0im

create(N::Int)

Bosonic creation operator with Hilbert space cutoff N.

This operator acts on a fock state as $\hat{a}^\dagger \ket{n} = \sqrt{n+1} \ket{n+1}$.

Examples

julia> a_d = create(20)
Quantum Object:   type=Operator   dims=[20]   size=(20, 20)   ishermitian=false
20×20 SparseMatrixCSC{ComplexF64, Int64} with 19 stored entries:
⎡⠢⡀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠈⠢⡀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠈⠢⡀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠈⠢⡀⠀⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠈⠢⡀⎦

julia> fock(20, 4)' * a_d * fock(20, 3)
2.0 + 0.0im

displace(N::Int, α::Number)

Generate a displacement operator:

\[\hat{D}(\alpha)=\exp\left( \alpha \hat{a}^\dagger - \alpha^* \hat{a} \right),\]

where $\hat{a}$ is the bosonic annihilation operator, and $\alpha$ is the amount of displacement in optical phase space.

squeeze(N::Int, z::Number)

Generate a single-mode squeeze operator:

\[\hat{S}(z)=\exp\left( \frac{1}{2} (z^* \hat{a}^2 - z(\hat{a}^\dagger)^2) \right),\]

where $\hat{a}$ is the bosonic annihilation operator.

num(N::Int)

Bosonic number operator with Hilbert space cutoff N.

This operator is defined as $\hat{N}=\hat{a}^\dagger \hat{a}$, where $\hat{a}$ is the bosonic annihilation operator.

position(N::Int)

Position operator with Hilbert space cutoff N.

This operator is defined as $\hat{x}=\frac{1}{\sqrt{2}} (\hat{a}^\dagger + \hat{a})$, where $\hat{a}$ is the bosonic annihilation operator.

momentum(N::Int)

Momentum operator with Hilbert space cutoff N.

This operator is defined as $\hat{p}= \frac{i}{\sqrt{2}} (\hat{a}^\dagger - \hat{a})$, where $\hat{a}$ is the bosonic annihilation operator.

phase(N::Int, ϕ0::Real=0)

Single-mode Pegg-Barnett phase operator with Hilbert space cutoff $N$ and the reference phase $\phi_0$.

This operator is defined as

\[\hat{\phi} = \sum_{m=0}^{N-1} \phi_m |\phi_m\rangle \langle\phi_m|,\]

where

\[\phi_m = \phi_0 + \frac{2m\pi}{N},\]

and

\[|\phi_m\rangle = \frac{1}{\sqrt{N}} \sum_{n=0}^{N-1} \exp(i n \phi_m) |n\rangle.\]

Reference

• Michael Martin Nieto, QUANTUM PHASE AND QUANTUM PHASE OPERATORS: Some Physics and Some History, arXiv:hep-th/9304036, Equation (30-32).

fdestroy(N::Union{Int,Val}, j::Int)

Construct a fermionic destruction operator acting on the j-th site, where the fock space has totally N-sites:

Here, we use the Jordan-Wigner transformation, namely

\[\hat{d}_j = \hat{\sigma}_z^{\otimes j-1} \otimes \hat{\sigma}_{+} \otimes \hat{\mathbb{1}}^{\otimes N-j}\]

The site index j should satisfy: 1 ≤ j ≤ N.

Note that we put $\hat{\sigma}_{+} = \begin{pmatrix} 0 & 1 \\ 0 & 0 \end{pmatrix}$ here because we consider $|0\rangle = \begin{pmatrix} 1 \\ 0 \end{pmatrix}$ to be ground (vacant) state, and $|1\rangle = \begin{pmatrix} 0 \\ 1 \end{pmatrix}$ to be excited (occupied) state.

fcreate(N::Union{Int,Val}, j::Int)

Construct a fermionic creation operator acting on the j-th site, where the fock space has totally N-sites:

Here, we use the Jordan-Wigner transformation, namely

\[\hat{d}^\dagger_j = \hat{\sigma}_z^{\otimes j-1} \otimes \hat{\sigma}_{-} \otimes \hat{\mathbb{1}}^{\otimes N-j}\]

The site index j should satisfy: 1 ≤ j ≤ N.

Note that we put $\hat{\sigma}_{-} = \begin{pmatrix} 0 & 0 \\ 1 & 0 \end{pmatrix}$ here because we consider $|0\rangle = \begin{pmatrix} 1 \\ 0 \end{pmatrix}$ to be ground (vacant) state, and $|1\rangle = \begin{pmatrix} 0 \\ 1 \end{pmatrix}$ to be excited (occupied) state.

tunneling(N::Int, m::Int=1; sparse::Union{Bool,Val{<:Bool}}=Val(false))

Generate a tunneling operator defined as:

\[\sum_{n=0}^{N-m} | n \rangle\langle n+m | + | n+m \rangle\langle n |,\]

where $N$ is the number of basis states in the Hilbert space, and $m$ is the number of excitations in tunneling event.

If sparse=true, the operator is returned as a sparse matrix, otherwise a dense matrix is returned.

qft(dimensions)

Generates a discrete Fourier transform matrix $\hat{F}_N$ for Quantum Fourier Transform (QFT) with given argument dimensions.

The dimensions can be either the following types:

• dimensions::Int: Number of basis states in the Hilbert space.
• dimensions::Union{AbstractVector{Int},Tuple}: list of dimensions representing the each number of basis in the subsystems.

$N$ represents the total dimension, and therefore the matrix is defined as

\[\hat{F}_N = \frac{1}{\sqrt{N}}\begin{bmatrix} 1 & 1 & 1 & 1 & \cdots & 1\\ 1 & \omega & \omega^2 & \omega^3 & \cdots & \omega^{N-1}\\ 1 & \omega^2 & \omega^4 & \omega^6 & \cdots & \omega^{2(N-1)}\\ 1 & \omega^3 & \omega^6 & \omega^9 & \cdots & \omega^{3(N-1)}\\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & \omega^{N-1} & \omega^{2(N-1)} & \omega^{3(N-1)} & \cdots & \omega^{(N-1)(N-1)} \end{bmatrix},\]

where $\omega = \exp(\frac{2 \pi i}{N})$.

eye(N::Int; type=Operator, dims=nothing)

Identity operator $\hat{\mathbb{1}}$ with size N.

It is also possible to specify the list of Hilbert dimensions dims if different subsystems are present.

Note that type can only be either Operator or SuperOperator

projection(N::Int, i::Int, j::Int)

Generates the projection operator $\hat{O} = \dyad{i}{j}$ with Hilbert space dimension N.

commutator(A::QuantumObject, B::QuantumObject; anti::Bool=false)

Return the commutator (or anti-commutator) of the two QuantumObject:

• commutator (anti=false): $\hat{A}\hat{B}-\hat{B}\hat{A}$
• anticommutator (anti=true): $\hat{A}\hat{B}+\hat{B}\hat{A}$

Note that A and B must be Operator

spre(A::QuantumObject, Id_cache=I(size(A,1)))

Returns the SuperOperator form of A acting on the left of the density matrix operator: $\mathcal{O} \left(\hat{A}\right) \left[ \hat{\rho} \right] = \hat{A} \hat{\rho}$.

Since the density matrix is vectorized in OperatorKet form: $|\hat{\rho}\rangle\rangle$, this SuperOperator is always a matrix $\hat{\mathbb{1}} \otimes \hat{A}$, namely

\[\mathcal{O} \left(\hat{A}\right) \left[ \hat{\rho} \right] = \hat{\mathbb{1}} \otimes \hat{A} ~ |\hat{\rho}\rangle\rangle\]

(see the section in documentation: Superoperators and Vectorized Operators for more details)

The optional argument Id_cache can be used to pass a precomputed identity matrix. This can be useful when the same function is applied multiple times with a known Hilbert space dimension.

spost(B::QuantumObject, Id_cache=I(size(B,1)))

Returns the SuperOperator form of B acting on the right of the density matrix operator: $\mathcal{O} \left(\hat{B}\right) \left[ \hat{\rho} \right] = \hat{\rho} \hat{B}$.

Since the density matrix is vectorized in OperatorKet form: $|\hat{\rho}\rangle\rangle$, this SuperOperator is always a matrix $\hat{B}^T \otimes \hat{\mathbb{1}}$, namely

\[\mathcal{O} \left(\hat{B}\right) \left[ \hat{\rho} \right] = \hat{B}^T \otimes \hat{\mathbb{1}} ~ |\hat{\rho}\rangle\rangle\]

(see the section in documentation: Superoperators and Vectorized Operators for more details)

The optional argument Id_cache can be used to pass a precomputed identity matrix. This can be useful when the same function is applied multiple times with a known Hilbert space dimension.

sprepost(A::QuantumObject, B::QuantumObject)

Returns the SuperOperator form of A and B acting on the left and right of the density matrix operator, respectively: $\mathcal{O} \left( \hat{A}, \hat{B} \right) \left[ \hat{\rho} \right] = \hat{A} \hat{\rho} \hat{B}$.

Since the density matrix is vectorized in OperatorKet form: $|\hat{\rho}\rangle\rangle$, this SuperOperator is always a matrix $\hat{B}^T \otimes \hat{A}$, namely

\[\mathcal{O} \left(\hat{A}, \hat{B}\right) \left[ \hat{\rho} \right] = \hat{B}^T \otimes \hat{A} ~ |\hat{\rho}\rangle\rangle = \textrm{spre}(\hat{A}) * \textrm{spost}(\hat{B}) ~ |\hat{\rho}\rangle\rangle\]

(see the section in documentation: Superoperators and Vectorized Operators for more details)

See also spre and spost.

lindblad_dissipator(O::QuantumObject, Id_cache=I(size(O,1))

Returns the Lindblad SuperOperator defined as

\[\mathcal{D} \left( \hat{O} \right) \left[ \hat{\rho} \right] = \frac{1}{2} \left( 2 \hat{O} \hat{\rho} \hat{O}^\dagger - \hat{O}^\dagger \hat{O} \hat{\rho} - \hat{\rho} \hat{O}^\dagger \hat{O} \right)\]

The optional argument Id_cache can be used to pass a precomputed identity matrix. This can be useful when the same function is applied multiple times with a known Hilbert space dimension.

See also spre and spost.

Qobj(A::AbstractArray; type::QuantumObjectType, dims::Vector{Int})

Generate QuantumObject

Note that this functions is same as QuantumObject(A; type=type, dims=dims)

shape(A::QuantumObject)

Returns a tuple containing each dimensions of the array in the QuantumObject.

Note that this function is same as size(A)

isherm(A::QuantumObject)

Test whether the QuantumObject is Hermitian.

Note that this functions is same as ishermitian(A)

trans(A::QuantumObject)

Lazy matrix transpose of the QuantumObject.

Note that this function is same as transpose(A)

dag(A::QuantumObject)

Lazy adjoint (conjugate transposition) of the QuantumObject

Note that this function is same as adjoint(A)

matrix_element(i::QuantumObject, A::QuantumObject j::QuantumObject)

Compute the generalized dot product dot(i, A*j) between three QuantumObject: $\langle i | \hat{A} | j \rangle$

Note that this function is same as dot(i, A, j)

Supports the following inputs:

• A is in the type of Operator, with i and j are both Ket.
• A is in the type of SuperOperator, with i and j are both OperatorKet

unit(A::QuantumObject, p::Real)

Return normalized QuantumObject so that its p-norm equals to unity, i.e. norm(A, p) == 1.

Support for the following types of QuantumObject:

• If A is Ket or Bra, default p = 2
• If A is Operator, default p = 1

Note that this function is same as normalize(A, p)

Also, see norm about its definition for different types of QuantumObject.

sqrtm(A::QuantumObject)

Matrix square root of Operator type of QuantumObject

Note that for other types of QuantumObject use sprt(A) instead.

logm(A::QuantumObject)

Matrix logarithm of QuantumObject

Note that this function is same as log(A) and only supports for Operator and SuperOperator

expm(A::QuantumObject)

Matrix exponential of QuantumObject

Note that this function is same as exp(A) and only supports for Operator and SuperOperator

sinm(A::QuantumObject)

Matrix sine of QuantumObject, defined as

$\sin \left( \hat{A} \right) = \frac{e^{i \hat{A}} - e^{-i \hat{A}}}{2 i}$

Note that this function is same as sin(A) and only supports for Operator and SuperOperator

cosm(A::QuantumObject)

Matrix cosine of QuantumObject, defined as

$\cos \left( \hat{A} \right) = \frac{e^{i \hat{A}} + e^{-i \hat{A}}}{2}$

Note that this function is same as cos(A) and only supports for Operator and SuperOperator

tensor(A::QuantumObject, B::QuantumObject, ...)

Returns the Kronecker product $\hat{A} \otimes \hat{B} \otimes \cdots$.

Note that this function is same as kron(A, B, ...)

Examples

julia> x = sigmax()
Quantum Object:   type=Operator   dims=[2]   size=(2, 2)   ishermitian=true
2×2 SparseMatrixCSC{ComplexF64, Int64} with 2 stored entries:
     ⋅      1.0+0.0im
 1.0+0.0im      ⋅

julia> x_list = fill(x, 3);

julia> tensor(x_list...)
Quantum Object:   type=Operator   dims=[2, 2, 2]   size=(8, 8)   ishermitian=true
8×8 SparseMatrixCSC{ComplexF64, Int64} with 8 stored entries:
     ⋅          ⋅          ⋅      …      ⋅          ⋅      1.0+0.0im
     ⋅          ⋅          ⋅             ⋅      1.0+0.0im      ⋅
     ⋅          ⋅          ⋅         1.0+0.0im      ⋅          ⋅
     ⋅          ⋅          ⋅             ⋅          ⋅          ⋅
     ⋅          ⋅          ⋅             ⋅          ⋅          ⋅
     ⋅          ⋅      1.0+0.0im  …      ⋅          ⋅          ⋅
     ⋅      1.0+0.0im      ⋅             ⋅          ⋅          ⋅
 1.0+0.0im      ⋅          ⋅             ⋅          ⋅          ⋅

⊗(A::QuantumObject, B::QuantumObject)

Returns the Kronecker product $\hat{A} \otimes \hat{B}$.

Note that this function is same as kron(A, B)

Examples

julia> a = destroy(20)
Quantum Object:   type=Operator   dims=[20]   size=(20, 20)   ishermitian=false
20×20 SparseMatrixCSC{ComplexF64, Int64} with 19 stored entries:
⠈⠢⡀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠈⠢⡀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠈⠢⡀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠈⠢⡀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠈⠢

julia> a ⊗ a
Quantum Object:   type=Operator   dims=[20, 20]   size=(400, 400)   ishermitian=false
400×400 SparseMatrixCSC{ComplexF64, Int64} with 361 stored entries:
⠀⠀⠘⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀⠀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⢦⡀⠀
⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠙⠦

qeye(N::Int; type=Operator, dims=nothing)

Identity operator $\hat{\mathbb{1}}$ with size N.

It is also possible to specify the list of Hilbert dimensions dims if different subsystems are present.

Note that this function is same as eye(N, type=type, dims=dims), and type can only be either Operator or SuperOperator

struct TimeEvolutionSol

A structure storing the results and some information from solving time evolution.

Fields (Attributes)

• times::AbstractVector: The time list of the evolution.
• states::Vector{QuantumObject}: The list of result states.
• expect::Matrix: The expectation values corresponding to each time point in times.
• retcode: The return code from the solver.
• alg: The algorithm which is used during the solving process.
• abstol::Real: The absolute tolerance which is used during the solving process.
• reltol::Real: The relative tolerance which is used during the solving process.

struct TimeEvolutionMCSol

A structure storing the results and some information from solving quantum trajectories of the Monte Carlo wave function time evolution.

Fields (Attributes)

• ntraj::Int: Number of trajectories
• times::AbstractVector: The time list of the evolution.
• states::Vector{Vector{QuantumObject}}: The list of result states in each trajectory.
• expect::Matrix: The expectation values (averaging all trajectories) corresponding to each time point in times.
• expect_all::Array: The expectation values corresponding to each trajectory and each time point in times
• jump_times::Vector{Vector{Real}}: The time records of every quantum jump occurred in each trajectory.
• jump_which::Vector{Vector{Int}}: The indices of the jump operators in c_ops that describe the corresponding quantum jumps occurred in each trajectory.
• converged::Bool: Whether the solution is converged or not.
• alg: The algorithm which is used during the solving process.
• abstol::Real: The absolute tolerance which is used during the solving process.
• reltol::Real: The relative tolerance which is used during the solving process.

struct TimeEvolutionSSESol

A structure storing the results and some information from solving trajectories of the Stochastic Shrodinger equation time evolution.

Fields (Attributes)

• ntraj::Int: Number of trajectories
• times::AbstractVector: The time list of the evolution.
• states::Vector{Vector{QuantumObject}}: The list of result states in each trajectory.
• expect::Matrix: The expectation values (averaging all trajectories) corresponding to each time point in times.
• expect_all::Array: The expectation values corresponding to each trajectory and each time point in times
• converged::Bool: Whether the solution is converged or not.
• alg: The algorithm which is used during the solving process.
• abstol::Real: The absolute tolerance which is used during the solving process.
• reltol::Real: The relative tolerance which is used during the solving process.

sesolveProblem(H::QuantumObject,
    ψ0::QuantumObject,
    tlist::AbstractVector;
    alg::OrdinaryDiffEqAlgorithm=Tsit5()
    e_ops::Union{Nothing,AbstractVector} = nothing,
    H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
    params::NamedTuple=NamedTuple(),
    progress_bar::Union{Val,Bool}=Val(true),
    kwargs...)

Generates the ODEProblem for the Schrödinger time evolution of a quantum system:

\[\frac{\partial}{\partial t} |\psi(t)\rangle = -i \hat{H} |\psi(t)\rangle\]

Arguments

• H::QuantumObject: The Hamiltonian of the system $\hat{H}$.
• ψ0::QuantumObject: The initial state of the system $|\psi(0)\rangle$.
• tlist::AbstractVector: The time list of the evolution.
• alg::OrdinaryDiffEqAlgorithm: The algorithm used for the time evolution.
• e_ops::Union{Nothing,AbstractVector}: The list of operators to be evaluated during the evolution.
• H_t::Union{Nothing,Function,TimeDependentOperatorSum}: The time-dependent Hamiltonian of the system. If nothing, the Hamiltonian is time-independent.
• params::NamedTuple: The parameters of the system.
• progress_bar::Union{Val,Bool}: Whether to show the progress bar. Using non-Val types might lead to type instabilities.
• kwargs...: The keyword arguments passed to the ODEProblem constructor.

Notes

• The states will be saved depend on the keyword argument saveat in kwargs.
• If e_ops is specified, the default value of saveat=[tlist[end]] (only save the final state), otherwise, saveat=tlist (saving the states corresponding to tlist). You can also specify e_ops and saveat separately.
• The default tolerances in kwargs are given as reltol=1e-6 and abstol=1e-8.
• For more details about alg please refer to DifferentialEquations.jl (ODE Solvers)
• For more details about kwargs please refer to DifferentialEquations.jl (Keyword Arguments)

Returns

• prob: The ODEProblem for the Schrödinger time evolution of the system.

mesolveProblem(H::QuantumObject,
    ψ0::QuantumObject,
    tlist::AbstractVector, 
    c_ops::Union{Nothing,AbstractVector}=nothing;
    alg::OrdinaryDiffEqAlgorithm=Tsit5(),
    e_ops::Union{Nothing,AbstractVector}=nothing,
    H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
    params::NamedTuple=NamedTuple(),
    progress_bar::Union{Val,Bool}=Val(true),
    kwargs...)

Generates the ODEProblem for the master equation time evolution of an open quantum system:

\[\frac{\partial \hat{\rho}(t)}{\partial t} = -i[\hat{H}, \hat{\rho}(t)] + \sum_n \mathcal{D}(\hat{C}_n) [\hat{\rho}(t)]\]

where

\[\mathcal{D}(\hat{C}_n) [\hat{\rho}(t)] = \hat{C}_n \hat{\rho}(t) \hat{C}_n^\dagger - \frac{1}{2} \hat{C}_n^\dagger \hat{C}_n \hat{\rho}(t) - \frac{1}{2} \hat{\rho}(t) \hat{C}_n^\dagger \hat{C}_n\]

Arguments

• H::QuantumObject: The Hamiltonian $\hat{H}$ or the Liouvillian of the system.
• ψ0::QuantumObject: The initial state of the system.
• tlist::AbstractVector: The time list of the evolution.
• c_ops::Union{Nothing,AbstractVector}=nothing: The list of the collapse operators $\{\hat{C}_n\}_n$.
• alg::OrdinaryDiffEqAlgorithm=Tsit5(): The algorithm used for the time evolution.
• e_ops::Union{Nothing,AbstractVector}=nothing: The list of the operators for which the expectation values are calculated.
• H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing: The time-dependent Hamiltonian or Liouvillian.
• params::NamedTuple=NamedTuple(): The parameters of the time evolution.
• progress_bar::Union{Val,Bool}=Val(true): Whether to show the progress bar. Using non-Val types might lead to type instabilities.
• kwargs...: The keyword arguments for the ODEProblem.

Notes

• The states will be saved depend on the keyword argument saveat in kwargs.
• If e_ops is specified, the default value of saveat=[tlist[end]] (only save the final state), otherwise, saveat=tlist (saving the states corresponding to tlist). You can also specify e_ops and saveat separately.
• The default tolerances in kwargs are given as reltol=1e-6 and abstol=1e-8.
• For more details about alg please refer to DifferentialEquations.jl (ODE Solvers)
• For more details about kwargs please refer to DifferentialEquations.jl (Keyword Arguments)

Returns

• prob::ODEProblem: The ODEProblem for the master equation time evolution.

mcsolveProblem(H::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
    ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
    tlist::AbstractVector,
    c_ops::Union{Nothing,AbstractVector}=nothing;
    alg::OrdinaryDiffEqAlgorithm=Tsit5(),
    e_ops::Union{Nothing,AbstractVector}=nothing,
    H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
    params::NamedTuple=NamedTuple(),
    jump_callback::TJC=ContinuousLindbladJumpCallback(),
    kwargs...)

Generates the ODEProblem for a single trajectory of the Monte Carlo wave function time evolution of an open quantum system.

Given a system Hamiltonian $\hat{H}$ and a list of collapse (jump) operators $\{\hat{C}_n\}_n$, the evolution of the state $|\psi(t)\rangle$ is governed by the Schrodinger equation:

\[\frac{\partial}{\partial t} |\psi(t)\rangle= -i \hat{H}_{\textrm{eff}} |\psi(t)\rangle\]

with a non-Hermitian effective Hamiltonian:

\[\hat{H}_{\textrm{eff}} = \hat{H} - \frac{i}{2} \sum_n \hat{C}_n^\dagger \hat{C}_n.\]

To the first-order of the wave function in a small time $\delta t$, the strictly negative non-Hermitian portion in $\hat{H}_{\textrm{eff}}$ gives rise to a reduction in the norm of the wave function, namely

\[\langle \psi(t+\delta t) | \psi(t+\delta t) \rangle = 1 - \delta p,\]

where

\[\delta p = \delta t \sum_n \langle \psi(t) | \hat{C}_n^\dagger \hat{C}_n | \psi(t) \rangle\]

is the corresponding quantum jump probability.

If the environmental measurements register a quantum jump, the wave function undergoes a jump into a state defined by projecting $|\psi(t)\rangle$ using the collapse operator $\hat{C}_n$ corresponding to the measurement, namely

\[| \psi(t+\delta t) \rangle = \frac{\hat{C}_n |\psi(t)\rangle}{ \sqrt{\langle \psi(t) | \hat{C}_n^\dagger \hat{C}_n | \psi(t) \rangle} }\]

Arguments

• H::QuantumObject: Hamiltonian of the system $\hat{H}$.
• ψ0::QuantumObject: Initial state of the system $|\psi(0)\rangle$.
• tlist::AbstractVector: List of times at which to save the state of the system.
• c_ops::Union{Nothing,AbstractVector}: List of collapse operators $\{\hat{C}_n\}_n$.
• alg::OrdinaryDiffEqAlgorithm: Algorithm to use for the time evolution.
• e_ops::Union{Nothing,AbstractVector}: List of operators for which to calculate expectation values.
• H_t::Union{Nothing,Function,TimeDependentOperatorSum}: Time-dependent part of the Hamiltonian.
• params::NamedTuple: Dictionary of parameters to pass to the solver.
• seeds::Union{Nothing, Vector{Int}}: List of seeds for the random number generator. Length must be equal to the number of trajectories provided.
• jump_callback::LindbladJumpCallbackType: The Jump Callback type: Discrete or Continuous.
• kwargs...: Additional keyword arguments to pass to the solver.

Notes

• The states will be saved depend on the keyword argument saveat in kwargs.
• If e_ops is specified, the default value of saveat=[tlist[end]] (only save the final state), otherwise, saveat=tlist (saving the states corresponding to tlist). You can also specify e_ops and saveat separately.
• The default tolerances in kwargs are given as reltol=1e-6 and abstol=1e-8.
• For more details about alg please refer to DifferentialEquations.jl (ODE Solvers)
• For more details about kwargs please refer to DifferentialEquations.jl (Keyword Arguments)

Returns

• prob::ODEProblem: The ODEProblem for the Monte Carlo wave function time evolution.

mcsolveEnsembleProblem(H::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
    ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
    tlist::AbstractVector,
    c_ops::Union{Nothing,AbstractVector}=nothing;
    alg::OrdinaryDiffEqAlgorithm=Tsit5(),
    e_ops::Union{Nothing,AbstractVector}=nothing,
    H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
    params::NamedTuple=NamedTuple(),
    jump_callback::TJC=ContinuousLindbladJumpCallback(),
    prob_func::Function=_mcsolve_prob_func,
    output_func::Function=_mcsolve_output_func,
    kwargs...)

Generates the EnsembleProblem of ODEProblems for the ensemble of trajectories of the Monte Carlo wave function time evolution of an open quantum system.

Given a system Hamiltonian $\hat{H}$ and a list of collapse (jump) operators $\{\hat{C}_n\}_n$, the evolution of the state $|\psi(t)\rangle$ is governed by the Schrodinger equation:

\[\frac{\partial}{\partial t} |\psi(t)\rangle= -i \hat{H}_{\textrm{eff}} |\psi(t)\rangle\]

with a non-Hermitian effective Hamiltonian:

\[\hat{H}_{\textrm{eff}} = \hat{H} - \frac{i}{2} \sum_n \hat{C}_n^\dagger \hat{C}_n.\]

To the first-order of the wave function in a small time $\delta t$, the strictly negative non-Hermitian portion in $\hat{H}_{\textrm{eff}}$ gives rise to a reduction in the norm of the wave function, namely

\[\langle \psi(t+\delta t) | \psi(t+\delta t) \rangle = 1 - \delta p,\]

where

\[\delta p = \delta t \sum_n \langle \psi(t) | \hat{C}_n^\dagger \hat{C}_n | \psi(t) \rangle\]

is the corresponding quantum jump probability.

If the environmental measurements register a quantum jump, the wave function undergoes a jump into a state defined by projecting $|\psi(t)\rangle$ using the collapse operator $\hat{C}_n$ corresponding to the measurement, namely

\[| \psi(t+\delta t) \rangle = \frac{\hat{C}_n |\psi(t)\rangle}{ \sqrt{\langle \psi(t) | \hat{C}_n^\dagger \hat{C}_n | \psi(t) \rangle} }\]

Arguments

• H::QuantumObject: Hamiltonian of the system $\hat{H}$.
• ψ0::QuantumObject: Initial state of the system $|\psi(0)\rangle$.
• tlist::AbstractVector: List of times at which to save the state of the system.
• c_ops::Union{Nothing,AbstractVector}: List of collapse operators $\{\hat{C}_n\}_n$.
• alg::OrdinaryDiffEqAlgorithm: Algorithm to use for the time evolution.
• e_ops::Union{Nothing,AbstractVector}: List of operators for which to calculate expectation values.
• H_t::Union{Nothing,Function,TimeDependentOperatorSum}: Time-dependent part of the Hamiltonian.
• params::NamedTuple: Dictionary of parameters to pass to the solver.
• seeds::Union{Nothing, Vector{Int}}: List of seeds for the random number generator. Length must be equal to the number of trajectories provided.
• jump_callback::LindbladJumpCallbackType: The Jump Callback type: Discrete or Continuous.
• prob_func::Function: Function to use for generating the ODEProblem.
• output_func::Function: Function to use for generating the output of a single trajectory.
• kwargs...: Additional keyword arguments to pass to the solver.

Notes

• The states will be saved depend on the keyword argument saveat in kwargs.
• If e_ops is specified, the default value of saveat=[tlist[end]] (only save the final state), otherwise, saveat=tlist (saving the states corresponding to tlist). You can also specify e_ops and saveat separately.
• The default tolerances in kwargs are given as reltol=1e-6 and abstol=1e-8.
• For more details about alg please refer to DifferentialEquations.jl (ODE Solvers)
• For more details about kwargs please refer to DifferentialEquations.jl (Keyword Arguments)

Returns

• prob::EnsembleProblem with ODEProblem: The Ensemble ODEProblem for the Monte Carlo wave function time evolution.

ssesolveProblem(H::QuantumObject,
    ψ0::QuantumObject,
    tlist::AbstractVector;
    sc_ops::Union{Nothing,AbstractVector}=nothing;
    alg::StochasticDiffEqAlgorithm=SRA1()
    e_ops::Union{Nothing,AbstractVector} = nothing,
    H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
    params::NamedTuple=NamedTuple(),
    progress_bar::Union{Val,Bool}=Val(true),
    kwargs...)

Generates the SDEProblem for the Stochastic Schrödinger time evolution of a quantum system. This is defined by the following stochastic differential equation:

\[ d|\psi(t)\rangle = -i K |\psi(t)\rangle dt + \sum_n M_n |\psi(t)\rangle dW_n(t) ``` where \]

math K = \hat{H} + i \sumn \left(\frac{ej} Cn - \frac{1}{2} \sum{j} Cn^\dagger Cn - \frac{ej^2}{8}\right), math Mn = Cn - \frac{en}{2}, andmath en = \langle Cn + C_n^\dagger \rangle. ```

Above, C_n is the n-th collapse operator and dW_j(t) is the real Wiener increment associated to C_n.

Arguments

• H::QuantumObject: The Hamiltonian of the system $\hat{H}$.
• ψ0::QuantumObject: The initial state of the system $|\psi(0)\rangle$.
• tlist::AbstractVector: The time list of the evolution.
• sc_ops::Union{Nothing,AbstractVector}=nothing: List of stochastic collapse operators $\{\hat{C}_n\}_n$.
• alg::StochasticDiffEqAlgorithm: The algorithm used for the time evolution.
• e_ops::Union{Nothing,AbstractVector}=nothing: The list of operators to be evaluated during the evolution.
• H_t::Union{Nothing,Function,TimeDependentOperatorSum}: The time-dependent Hamiltonian of the system. If nothing, the Hamiltonian is time-independent.
• params::NamedTuple: The parameters of the system.
• progress_bar::Union{Val,Bool}: Whether to show the progress bar. Using non-Val types might lead to type instabilities.
• kwargs...: The keyword arguments passed to the SDEProblem constructor.

Notes

• The states will be saved depend on the keyword argument saveat in kwargs.
• If e_ops is specified, the default value of saveat=[tlist[end]] (only save the final state), otherwise, saveat=tlist (saving the states corresponding to tlist). You can also specify e_ops and saveat separately.
• The default tolerances in kwargs are given as reltol=1e-2 and abstol=1e-2.
• For more details about alg please refer to DifferentialEquations.jl (SDE Solvers)
• For more details about kwargs please refer to DifferentialEquations.jl (Keyword Arguments)

Returns

• prob: The SDEProblem for the Stochastic Schrödinger time evolution of the system.

ssesolveEnsembleProblem(H::QuantumObject,
    ψ0::QuantumObject,
    tlist::AbstractVector;
    sc_ops::Union{Nothing,AbstractVector} = nothing;
    alg::StochasticDiffEqAlgorithm=SRA1()
    e_ops::Union{Nothing,AbstractVector} = nothing,
    H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
    params::NamedTuple=NamedTuple(),
    prob_func::Function=_mcsolve_prob_func,
    output_func::Function=_mcsolve_output_func,
    kwargs...)

Generates the SDE EnsembleProblem for the Stochastic Schrödinger time evolution of a quantum system. This is defined by the following stochastic differential equation:

\[ d|\psi(t)\rangle = -i K |\psi(t)\rangle dt + \sum_n M_n |\psi(t)\rangle dW_n(t) ``` where \]

math K = \hat{H} + i \sumn \left(\frac{ej} Cn - \frac{1}{2} \sum{j} Cn^\dagger Cn - \frac{ej^2}{8}\right), math Mn = Cn - \frac{en}{2}, andmath en = \langle Cn + C_n^\dagger \rangle. ```

Above, C_n is the n-th collapse operator and dW_j(t) is the real Wiener increment associated to C_n.

Arguments

• H::QuantumObject: The Hamiltonian of the system $\hat{H}$.
• ψ0::QuantumObject: The initial state of the system $|\psi(0)\rangle$.
• tlist::AbstractVector: The time list of the evolution.
• sc_ops::Union{Nothing,AbstractVector}=nothing: List of stochastic collapse operators $\{\hat{C}_n\}_n$.
• alg::StochasticDiffEqAlgorithm: The algorithm used for the time evolution.
• e_ops::Union{Nothing,AbstractVector}=nothing: The list of operators to be evaluated during the evolution.
• H_t::Union{Nothing,Function,TimeDependentOperatorSum}: The time-dependent Hamiltonian of the system. If nothing, the Hamiltonian is time-independent.
• params::NamedTuple: The parameters of the system.
• prob_func::Function: Function to use for generating the SDEProblem.
• output_func::Function: Function to use for generating the output of a single trajectory.
• kwargs...: The keyword arguments passed to the SDEProblem constructor.

Notes

• The states will be saved depend on the keyword argument saveat in kwargs.
• If e_ops is specified, the default value of saveat=[tlist[end]] (only save the final state), otherwise, saveat=tlist (saving the states corresponding to tlist). You can also specify e_ops and saveat separately.
• The default tolerances in kwargs are given as reltol=1e-2 and abstol=1e-2.
• For more details about alg please refer to DifferentialEquations.jl (SDE Solvers)
• For more details about kwargs please refer to DifferentialEquations.jl (Keyword Arguments)

Returns

• prob::EnsembleProblem with SDEProblem: The Ensemble SDEProblem for the Stochastic Shrödinger time evolution.

lr_mesolveProblem(H, z, B, tlist, c_ops; e_ops=(), f_ops=(), opt=LRMesolveOptions(), kwargs...) where T
Formulates the ODEproblem for the low-rank time evolution of the system. The function is called by lr_mesolve.

Parameters
----------
H : QuantumObject
    The Hamiltonian of the system.
z : AbstractMatrix{T}
    The initial z matrix.
B : AbstractMatrix{T}
    The initial B matrix.
tlist : AbstractVector{T}
    The time steps at which the expectation values and function values are calculated.
c_ops : AbstractVector{QuantumObject}
    The jump operators of the system.
e_ops : Tuple{QuantumObject}
    The operators whose expectation values are calculated.
f_ops : Tuple{Function}
    The functions whose values are calculated.
opt : LRMesolveOptions
    The options of the problem.
kwargs : NamedTuple
    Additional keyword arguments for the ODEProblem.

sesolve(H::QuantumObject,
    ψ0::QuantumObject,
    tlist::AbstractVector;
    alg::OrdinaryDiffEqAlgorithm=Tsit5(),
    e_ops::Union{Nothing,AbstractVector} = nothing,
    H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
    params::NamedTuple=NamedTuple(),
    progress_bar::Union{Val,Bool}=Val(true),
    kwargs...)

Time evolution of a closed quantum system using the Schrödinger equation:

\[\frac{\partial}{\partial t} |\psi(t)\rangle = -i \hat{H} |\psi(t)\rangle\]

Arguments

• H::QuantumObject: The Hamiltonian of the system $\hat{H}$.
• ψ0::QuantumObject: The initial state of the system $|\psi(0)\rangle$.
• tlist::AbstractVector: List of times at which to save the state of the system.
• alg::OrdinaryDiffEqAlgorithm: Algorithm to use for the time evolution.
• e_ops::Union{Nothing,AbstractVector}: List of operators for which to calculate expectation values.
• H_t::Union{Nothing,Function,TimeDependentOperatorSum}: Time-dependent part of the Hamiltonian.
• params::NamedTuple: Dictionary of parameters to pass to the solver.
• progress_bar::Union{Val,Bool}: Whether to show the progress bar. Using non-Val types might lead to type instabilities.
• kwargs...: Additional keyword arguments to pass to the solver.

Notes

• The states will be saved depend on the keyword argument saveat in kwargs.
• If e_ops is specified, the default value of saveat=[tlist[end]] (only save the final state), otherwise, saveat=tlist (saving the states corresponding to tlist). You can also specify e_ops and saveat separately.
• The default tolerances in kwargs are given as reltol=1e-6 and abstol=1e-8.
• For more details about alg please refer to DifferentialEquations.jl (ODE Solvers)
• For more details about kwargs please refer to DifferentialEquations.jl (Keyword Arguments)

Returns

• sol::TimeEvolutionSol: The solution of the time evolution. See also TimeEvolutionSol

mesolve(H::QuantumObject,
    ψ0::QuantumObject,
    tlist::AbstractVector, 
    c_ops::Union{Nothing,AbstractVector}=nothing;
    alg::OrdinaryDiffEqAlgorithm=Tsit5(),
    e_ops::Union{Nothing,AbstractVector}=nothing,
    H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
    params::NamedTuple=NamedTuple(),
    progress_bar::Union{Val,Bool}=Val(true),
    kwargs...)

Time evolution of an open quantum system using Lindblad master equation:

\[\frac{\partial \hat{\rho}(t)}{\partial t} = -i[\hat{H}, \hat{\rho}(t)] + \sum_n \mathcal{D}(\hat{C}_n) [\hat{\rho}(t)]\]

where

\[\mathcal{D}(\hat{C}_n) [\hat{\rho}(t)] = \hat{C}_n \hat{\rho}(t) \hat{C}_n^\dagger - \frac{1}{2} \hat{C}_n^\dagger \hat{C}_n \hat{\rho}(t) - \frac{1}{2} \hat{\rho}(t) \hat{C}_n^\dagger \hat{C}_n\]

Arguments

• H::QuantumObject: The Hamiltonian $\hat{H}$ or the Liouvillian of the system.
• ψ0::QuantumObject: The initial state of the system.
• tlist::AbstractVector: The time list of the evolution.
• c_ops::Union{Nothing,AbstractVector}=nothing: The list of the collapse operators $\{\hat{C}_n\}_n$.
• alg::OrdinaryDiffEqAlgorithm: Algorithm to use for the time evolution.
• e_ops::Union{Nothing,AbstractVector}: List of operators for which to calculate expectation values.
• H_t::Union{Nothing,Function,TimeDependentOperatorSum}: Time-dependent part of the Hamiltonian.
• params::NamedTuple: Named Tuple of parameters to pass to the solver.
• progress_bar::Union{Val,Bool}: Whether to show the progress bar. Using non-Val types might lead to type instabilities.
• kwargs...: Additional keyword arguments to pass to the solver.

Notes

• The states will be saved depend on the keyword argument saveat in kwargs.
• If e_ops is specified, the default value of saveat=[tlist[end]] (only save the final state), otherwise, saveat=tlist (saving the states corresponding to tlist). You can also specify e_ops and saveat separately.
• The default tolerances in kwargs are given as reltol=1e-6 and abstol=1e-8.
• For more details about alg please refer to DifferentialEquations.jl (ODE Solvers)
• For more details about kwargs please refer to DifferentialEquations.jl (Keyword Arguments)

Returns

• sol::TimeEvolutionSol: The solution of the time evolution. See also TimeEvolutionSol

mcsolve(H::QuantumObject{<:AbstractArray{T1},OperatorQuantumObject},
    ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
    tlist::AbstractVector,
    c_ops::Union{Nothing,AbstractVector}=nothing;
    alg::OrdinaryDiffEqAlgorithm=Tsit5(),
    e_ops::Union{Nothing,AbstractVector}=nothing,
    H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
    params::NamedTuple=NamedTuple(),
    ntraj::Int=1,
    ensemble_method=EnsembleThreads(),
    jump_callback::TJC=ContinuousLindbladJumpCallback(),
    kwargs...)

Time evolution of an open quantum system using quantum trajectories.

Given a system Hamiltonian $\hat{H}$ and a list of collapse (jump) operators $\{\hat{C}_n\}_n$, the evolution of the state $|\psi(t)\rangle$ is governed by the Schrodinger equation:

\[\frac{\partial}{\partial t} |\psi(t)\rangle= -i \hat{H}_{\textrm{eff}} |\psi(t)\rangle\]

with a non-Hermitian effective Hamiltonian:

\[\hat{H}_{\textrm{eff}} = \hat{H} - \frac{i}{2} \sum_n \hat{C}_n^\dagger \hat{C}_n.\]

To the first-order of the wave function in a small time $\delta t$, the strictly negative non-Hermitian portion in $\hat{H}_{\textrm{eff}}$ gives rise to a reduction in the norm of the wave function, namely

\[\langle \psi(t+\delta t) | \psi(t+\delta t) \rangle = 1 - \delta p,\]

where

\[\delta p = \delta t \sum_n \langle \psi(t) | \hat{C}_n^\dagger \hat{C}_n | \psi(t) \rangle\]

is the corresponding quantum jump probability.

If the environmental measurements register a quantum jump, the wave function undergoes a jump into a state defined by projecting $|\psi(t)\rangle$ using the collapse operator $\hat{C}_n$ corresponding to the measurement, namely

\[| \psi(t+\delta t) \rangle = \frac{\hat{C}_n |\psi(t)\rangle}{ \sqrt{\langle \psi(t) | \hat{C}_n^\dagger \hat{C}_n | \psi(t) \rangle} }\]

Arguments

• H::QuantumObject: Hamiltonian of the system $\hat{H}$.
• ψ0::QuantumObject: Initial state of the system $|\psi(0)\rangle$.
• tlist::AbstractVector: List of times at which to save the state of the system.
• c_ops::Union{Nothing,AbstractVector}: List of collapse operators $\{\hat{C}_n\}_n$.
• alg::OrdinaryDiffEqAlgorithm: Algorithm to use for the time evolution.
• e_ops::Union{Nothing,AbstractVector}: List of operators for which to calculate expectation values.
• H_t::Union{Nothing,Function,TimeDependentOperatorSum}: Time-dependent part of the Hamiltonian.
• params::NamedTuple: Dictionary of parameters to pass to the solver.
• seeds::Union{Nothing, Vector{Int}}: List of seeds for the random number generator. Length must be equal to the number of trajectories provided.
• ntraj::Int: Number of trajectories to use.
• ensemble_method: Ensemble method to use.
• jump_callback::LindbladJumpCallbackType: The Jump Callback type: Discrete or Continuous.
• prob_func::Function: Function to use for generating the ODEProblem.
• output_func::Function: Function to use for generating the output of a single trajectory.
• kwargs...: Additional keyword arguments to pass to the solver.

Notes

• ensemble_method can be one of EnsembleThreads(), EnsembleSerial(), EnsembleDistributed()
• The states will be saved depend on the keyword argument saveat in kwargs.
• If e_ops is specified, the default value of saveat=[tlist[end]] (only save the final state), otherwise, saveat=tlist (saving the states corresponding to tlist). You can also specify e_ops and saveat separately.
• The default tolerances in kwargs are given as reltol=1e-6 and abstol=1e-8.
• For more details about alg please refer to DifferentialEquations.jl (ODE Solvers)
• For more details about kwargs please refer to DifferentialEquations.jl (Keyword Arguments)

Returns

• sol::TimeEvolutionMCSol: The solution of the time evolution. See also TimeEvolutionMCSol

ssesolve(H::QuantumObject,
    ψ0::QuantumObject{<:AbstractArray{T2},KetQuantumObject},
    tlist::AbstractVector,
    sc_ops::Union{Nothing, AbstractVector}=nothing;
    alg::StochasticDiffEqAlgorithm=SRA1(),
    e_ops::Union{Nothing,AbstractVector}=nothing,
    H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
    params::NamedTuple=NamedTuple(),
    ntraj::Int=1,
    ensemble_method=EnsembleThreads(),
    prob_func::Function=_mcsolve_prob_func,
    output_func::Function=_mcsolve_output_func,
    kwargs...)

Stochastic Schrödinger equation evolution of a quantum system given the system Hamiltonian $\hat{H}$ and a list of stochadtic collapse (jump) operators $\{\hat{C}_n\}_n$. The stochastic evolution of the state $|\psi(t)\rangle$ is defined by:

\[ d|\psi(t)\rangle = -i K |\psi(t)\rangle dt + \sum_n M_n |\psi(t)\rangle dW_n(t) ``` where \]

math K = \hat{H} + i \sumn \left(\frac{ej} Cn - \frac{1}{2} \sum{j} Cn^\dagger Cn - \frac{ej^2}{8}\right), math Mn = Cn - \frac{en}{2}, andmath en = \langle Cn + C_n^\dagger \rangle. ```

Above, C_n is the n-th collapse operator and dW_j(t) is the real Wiener increment associated to C_n.

Arguments

• H::QuantumObject: Hamiltonian of the system $\hat{H}$.
• ψ0::QuantumObject: Initial state of the system $|\psi(0)\rangle$.
• tlist::AbstractVector: List of times at which to save the state of the system.
• sc_ops::Union{Nothing,AbstractVector}=nothing: List of stochastic collapse operators $\{\hat{C}_n\}_n$.
• alg::StochasticDiffEqAlgorithm: Algorithm to use for the time evolution.
• e_ops::Union{Nothing,AbstractVector}: List of operators for which to calculate expectation values.
• H_t::Union{Nothing,Function,TimeDependentOperatorSum}: Time-dependent part of the Hamiltonian.
• params::NamedTuple: Dictionary of parameters to pass to the solver.
• seeds::Union{Nothing, Vector{Int}}: List of seeds for the random number generator. Length must be equal to the number of trajectories provided.
• ntraj::Int: Number of trajectories to use.
• ensemble_method: Ensemble method to use.
• prob_func::Function: Function to use for generating the SDEProblem.
• output_func::Function: Function to use for generating the output of a single trajectory.
• kwargs...: Additional keyword arguments to pass to the solver.

Notes

• ensemble_method can be one of EnsembleThreads(), EnsembleSerial(), EnsembleDistributed()
• The states will be saved depend on the keyword argument saveat in kwargs.
• If e_ops is specified, the default value of saveat=[tlist[end]] (only save the final state), otherwise, saveat=tlist (saving the states corresponding to tlist). You can also specify e_ops and saveat separately.
• The default tolerances in kwargs are given as reltol=1e-2 and abstol=1e-2.
• For more details about alg please refer to DifferentialEquations.jl (SDE Solvers)
• For more details about kwargs please refer to DifferentialEquations.jl (Keyword Arguments)

Returns

• sol::TimeEvolutionSSESol: The solution of the time evolution. See also TimeEvolutionSSESol

dfd_mesolve(H::Function, ψ0::QuantumObject,
    t_l::AbstractVector, c_ops::Function, maxdims::AbstractVector,
    dfd_params::NamedTuple=NamedTuple();
    alg::OrdinaryDiffEqAlgorithm=Tsit5(),
    e_ops::Function=(dim_list) -> Vector{Vector{T1}}([]),
    H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
    params::NamedTuple=NamedTuple(),
    tol_list::Vector{<:Number}=fill(1e-8, length(maxdims)),
    kwargs...)

Time evolution of an open quantum system using master equation, dynamically changing the dimension of the Hilbert subspaces.

Notes

• For more details about alg please refer to DifferentialEquations.jl (ODE Solvers)
• For more details about kwargs please refer to DifferentialEquations.jl (Keyword Arguments)

dsf_mesolve(H::Function,
    ψ0::QuantumObject,
    t_l::AbstractVector, c_ops::Function,
    op_list::Vector{TOl},
    α0_l::Vector{<:Number}=zeros(length(op_list)),
    dsf_params::NamedTuple=NamedTuple();
    alg::OrdinaryDiffEqAlgorithm=Tsit5(),
    e_ops::Function=(op_list,p) -> Vector{TOl}([]),
    H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
    params::NamedTuple=NamedTuple(),
    δα_list::Vector{<:Number}=fill(0.2, length(op_list)),
    krylov_dim::Int=max(6, min(10, cld(length(ket2dm(ψ0).data), 4))),
    kwargs...)

Time evolution of an open quantum system using master equation and the Dynamical Shifted Fock algorithm.

Notes

• For more details about alg please refer to DifferentialEquations.jl (ODE Solvers)
• For more details about kwargs please refer to DifferentialEquations.jl (Keyword Arguments)

dsf_mcsolve(H::Function,
    ψ0::QuantumObject,
    t_l::AbstractVector, c_ops::Function,
    op_list::Vector{TOl},
    α0_l::Vector{<:Number}=zeros(length(op_list)),
    dsf_params::NamedTuple=NamedTuple();
    alg::OrdinaryDiffEqAlgorithm=Tsit5(),
    e_ops::Function=(op_list,p) -> Vector{TOl}([]),
    H_t::Union{Nothing,Function,TimeDependentOperatorSum}=nothing,
    params::NamedTuple=NamedTuple(),
    δα_list::Vector{<:Real}=fill(0.2, length(op_list)),
    ntraj::Int=1,
    ensemble_method=EnsembleThreads(),
    jump_callback::LindbladJumpCallbackType=ContinuousLindbladJumpCallback(),
    krylov_dim::Int=max(6, min(10, cld(length(ket2dm(ψ0).data), 4))),
    kwargs...)

Time evolution of a quantum system using the Monte Carlo wave function method and the Dynamical Shifted Fock algorithm.

Notes

• For more details about alg please refer to DifferentialEquations.jl (ODE Solvers)
• For more details about kwargs please refer to DifferentialEquations.jl (Keyword Arguments)

lr_mesolve(prob::ODEProblem; kwargs...)
Solves the ODEProblem formulated by lr_mesolveProblem. The function is called by lr_mesolve.

Parameters
----------
prob : ODEProblem
    The ODEProblem formulated by lr_mesolveProblem.
kwargs : NamedTuple
    Additional keyword arguments for the ODEProblem.

liouvillian(H::QuantumObject, c_ops::Union{AbstractVector,Nothing}=nothing, Id_cache=I(prod(H.dims)))

Construct the Liouvillian SuperOperator for a system Hamiltonian $\hat{H}$ and a set of collapse operators $\{\hat{C}_n\}_n$:

\[\mathcal{L} [\cdot] = -i[\hat{H}, \cdot] + \sum_n \mathcal{D}(\hat{C}_n) [\cdot]\]

where

\[\mathcal{D}(\hat{C}_n) [\cdot] = \hat{C}_n [\cdot] \hat{C}_n^\dagger - \frac{1}{2} \hat{C}_n^\dagger \hat{C}_n [\cdot] - \frac{1}{2} [\cdot] \hat{C}_n^\dagger \hat{C}_n\]

The optional argument Id_cache can be used to pass a precomputed identity matrix. This can be useful when the same function is applied multiple times with a known Hilbert space dimension.

See also spre, spost, and lindblad_dissipator.

liouvillian_generalized(H::QuantumObject, fields::Vector, 
T_list::Vector; N_trunc::Int=size(H,1), tol::Float64=0.0, σ_filter::Union{Nothing, Real}=nothing)

Constructs the generalized Liouvillian for a system coupled to a bath of harmonic oscillators.

See, e.g., Settineri, Alessio, et al. "Dissipation and thermal noise in hybrid quantum systems in the ultrastrong-coupling regime." Physical Review A 98.5 (2018): 053834.

steadystate(
    H::QuantumObject,
    c_ops::Union{Nothing,AbstractVector} = nothing;
    solver::SteadyStateSolver = SteadyStateDirectSolver(),
    kwargs...
)

Solve the stationary state based on different solvers.

Parameters

• H::QuantumObject: The Hamiltonian or the Liouvillian of the system.
• c_ops::Union{Nothing,AbstractVector}=nothing: The list of the collapse operators.
• solver::SteadyStateSolver=SteadyStateDirectSolver(): see documentation Solving for Steady-State Solutions for different solvers.
• kwargs...: The keyword arguments for the solver.

steadystate(
    H::QuantumObject,
    ψ0::QuantumObject,
    tspan::Real = Inf,
    c_ops::Union{Nothing,AbstractVector} = nothing;
    solver::SteadyStateODESolver = SteadyStateODESolver(),
    reltol::Real = 1.0e-8,
    abstol::Real = 1.0e-10,
    kwargs...
)

Solve the stationary state based on time evolution (ordinary differential equations; OrdinaryDiffEq.jl) with a given initial state.

The termination condition of the stationary state $|\rho\rangle\rangle$ is that either the following condition is true:

\[\lVert\frac{\partial |\hat{\rho}\rangle\rangle}{\partial t}\rVert \leq \textrm{reltol} \times\lVert\frac{\partial |\hat{\rho}\rangle\rangle}{\partial t}+|\hat{\rho}\rangle\rangle\rVert\]

or

\[\lVert\frac{\partial |\hat{\rho}\rangle\rangle}{\partial t}\rVert \leq \textrm{abstol}\]

Parameters

• H::QuantumObject: The Hamiltonian or the Liouvillian of the system.
• ψ0::QuantumObject: The initial state of the system.
• tspan::Real=Inf: The final time step for the steady state problem.
• c_ops::Union{Nothing,AbstractVector}=nothing: The list of the collapse operators.
• solver::SteadyStateODESolver=SteadyStateODESolver(): see SteadyStateODESolver for more details.
• reltol::Real=1.0e-8: Relative tolerance in steady state terminate condition and solver adaptive timestepping.
• abstol::Real=1.0e-10: Absolute tolerance in steady state terminate condition and solver adaptive timestepping.
• kwargs...: The keyword arguments for the ODEProblem.

steadystate_floquet(
    H_0::QuantumObject{MT,OpType1},
    H_p::QuantumObject{<:AbstractArray,OpType2},
    H_m::QuantumObject{<:AbstractArray,OpType3},
    ωd::Number,
    c_ops::Union{Nothing,AbstractVector} = nothing;
    n_max::Integer = 2,
    tol::R = 1e-8,
    solver::FSolver = SSFloquetLinearSystem,
    kwargs...,
)

Calculates the steady state of a periodically driven system. Here H_0 is the Hamiltonian or the Liouvillian of the undriven system. Considering a monochromatic drive at frequency $\omega_d$, we divide it into two parts, H_p and H_m, where H_p oscillates as $e^{i \omega t}$ and H_m oscillates as $e^{-i \omega t}$. There are two solvers available for this function:

• SSFloquetLinearSystem: Solves the linear system of equations.
• SSFloquetEffectiveLiouvillian: Solves the effective Liouvillian.

For both cases, n_max is the number of Fourier components to consider, and tol is the tolerance for the solver.

In the case of SSFloquetLinearSystem, the full linear system is solved at once:

\[( \mathcal{L}_0 - i n \omega_d ) \hat{\rho}_n + \mathcal{L}_1 \hat{\rho}_{n-1} + \mathcal{L}_{-1} \hat{\rho}_{n+1} = 0\]

This is a tridiagonal linear system of the form

\[\mathbf{A} \cdot \mathbf{b} = 0\]

where

\[\mathbf{A} = \begin{pmatrix} \mathcal{L}_0 - i (-n_{\textrm{max}}) \omega_{\textrm{d}} & \mathcal{L}_{-1} & 0 & \cdots & 0 \\ \mathcal{L}_1 & \mathcal{L}_0 - i (-n_{\textrm{max}}+1) \omega_{\textrm{d}} & \mathcal{L}_{-1} & \cdots & 0 \\ 0 & \mathcal{L}_1 & \mathcal{L}_0 - i (-n_{\textrm{max}}+2) \omega_{\textrm{d}} & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & \mathcal{L}_0 - i n_{\textrm{max}} \omega_{\textrm{d}} \end{pmatrix}\]

and

\[\mathbf{b} = \begin{pmatrix} \hat{\rho}_{-n_{\textrm{max}}} \\ \hat{\rho}_{-n_{\textrm{max}}+1} \\ \vdots \\ \hat{\rho}_{0} \\ \vdots \\ \hat{\rho}_{n_{\textrm{max}}-1} \\ \hat{\rho}_{n_{\textrm{max}}} \end{pmatrix}\]

This will allow to simultaneously obtain all the $\hat{\rho}_n$.

In the case of SSFloquetEffectiveLiouvillian, instead, the effective Liouvillian is calculated using the matrix continued fraction method.

Arguments

• H_0::QuantumObject: The Hamiltonian or the Liouvillian of the undriven system.
• H_p::QuantumObject: The Hamiltonian or the Liouvillian of the part of the drive that oscillates as $e^{i \omega t}$.
• H_m::QuantumObject: The Hamiltonian or the Liouvillian of the part of the drive that oscillates as $e^{-i \omega t}$.
• ωd::Number: The frequency of the drive.
• c_ops::AbstractVector = QuantumObject: The optional collapse operators.
• n_max::Integer = 2: The number of Fourier components to consider.
• tol::R = 1e-8: The tolerance for the solver.
• solver::FSolver = SSFloquetLinearSystem: The solver to use.
• kwargs...: Additional keyword arguments to be passed to the solver.

SteadyStateDirectSolver()

A solver which solves steadystate by finding the inverse of Liouvillian SuperOperator using the standard method given in LinearAlgebra.

SteadyStateEigenSolver()

A solver which solves steadystate by finding the zero (or lowest) eigenvalue of Liouvillian SuperOperator using eigsolve.

SteadyStateLinearSolver(alg = KrylovJL_GMRES(), Pl = nothing, Pr = nothing)

A solver which solves steadystate by finding the inverse of Liouvillian SuperOperator using the algorithms given in LinearSolve.jl.

Parameters

• alg::SciMLLinearSolveAlgorithm=KrylovJL_GMRES(): algorithms given in LinearSolve.jl
• Pl::Union{Function,Nothing}=nothing: left preconditioner, see documentation Solving for Steady-State Solutions for more details.
• Pr::Union{Function,Nothing}=nothing: right preconditioner, see documentation Solving for Steady-State Solutions for more details.

SteadyStateODESolver(alg = Tsit5())

An ordinary differential equation (ODE) solver for solving steadystate.

It includes a field (attribute) SteadyStateODESolver.alg that specifies the solving algorithm. Default to Tsit5().

For more details about the solvers, please refer to OrdinaryDiffEq.jl

correlation_3op_2t(H::QuantumObject,
    ψ0::QuantumObject,
    t_l::AbstractVector,
    τ_l::AbstractVector,
    A::QuantumObject,
    B::QuantumObject,
    C::QuantumObject,
    c_ops::Union{Nothing,AbstractVector}=nothing;
    kwargs...)

Returns the two-times correlation function of three operators $\hat{A}$, $\hat{B}$ and $\hat{C}$: $\expval{\hat{A}(t) \hat{B}(t + \tau) \hat{C}(t)}$

for a given initial state $\ket{\psi_0}$.

correlation_2op_2t(H::QuantumObject,
    ψ0::QuantumObject,
    t_l::AbstractVector,
    τ_l::AbstractVector,
    A::QuantumObject,
    B::QuantumObject,
    c_ops::Union{Nothing,AbstractVector}=nothing;
    reverse::Bool=false,
    kwargs...)

Returns the two-times correlation function of two operators $\hat{A}$ and $\hat{B}$ at different times: $\expval{\hat{A}(t + \tau) \hat{B}(t)}$.

When reverse=true, the correlation function is calculated as $\expval{\hat{A}(t) \hat{B}(t + \tau)}$.

correlation_2op_1t(H::QuantumObject,
    ψ0::QuantumObject,
    τ_l::AbstractVector,
    A::QuantumObject,
    B::QuantumObject,
    c_ops::Union{Nothing,AbstractVector}=nothing;
    reverse::Bool=false,
    kwargs...)

Returns the one-time correlation function of two operators $\hat{A}$ and $\hat{B}$ at different times $\expval{\hat{A}(\tau) \hat{B}(0)}$.

When reverse=true, the correlation function is calculated as $\expval{\hat{A}(0) \hat{B}(\tau)}$.

spectrum(H::QuantumObject,
    ω_list::AbstractVector,
    A::QuantumObject{<:AbstractArray{T2},OperatorQuantumObject},
    B::QuantumObject{<:AbstractArray{T3},OperatorQuantumObject},
    c_ops::Union{Nothing,AbstractVector}=nothing;
    solver::MySolver=ExponentialSeries(),
    kwargs...)

Returns the emission spectrum

\[S(\omega) = \int_{-\infty}^\infty \expval{\hat{A}(\tau) \hat{B}(0)} e^{-i \omega \tau} d \tau\]

entropy_vn(ρ::QuantumObject; base::Int=0, tol::Real=1e-15)

Calculates the Von Neumann entropy $S = - \Tr \left[ \hat{\rho} \log \left( \hat{\rho} \right) \right]$ where $\hat{\rho}$ is the density matrix of the system.

The base parameter specifies the base of the logarithm to use, and when using the default value 0, the natural logarithm is used. The tol parameter describes the absolute tolerance for detecting the zero-valued eigenvalues of the density matrix $\hat{\rho}$.

Examples

Pure state:

julia> ψ = fock(2,0)
Quantum Object:   type=Ket   dims=[2]   size=(2,)
2-element Vector{ComplexF64}:
 1.0 + 0.0im
 0.0 + 0.0im

julia> ρ = ket2dm(ψ)
Quantum Object:   type=Operator   dims=[2]   size=(2, 2)   ishermitian=true
2×2 Matrix{ComplexF64}:
 1.0+0.0im  0.0+0.0im
 0.0+0.0im  0.0+0.0im

julia> entropy_vn(ρ, base=2)
-0.0

Mixed state:

julia> ρ = maximally_mixed_dm(2)
Quantum Object:   type=Operator   dims=[2]   size=(2, 2)   ishermitian=true
2×2 Diagonal{ComplexF64, Vector{ComplexF64}}:
 0.5-0.0im      ⋅    
     ⋅      0.5-0.0im

julia> entropy_vn(ρ, base=2)
1.0

entanglement(QO::QuantumObject, sel::Union{Int,AbstractVector{Int},Tuple})

Calculates the entanglement by doing the partial trace of QO, selecting only the dimensions with the indices contained in the sel vector, and then using the Von Neumann entropy entropy_vn.

tracedist(ρ::QuantumObject, σ::QuantumObject)

Calculates the trace distance between two QuantumObject: $T(\hat{\rho}, \hat{\sigma}) = \frac{1}{2} \lVert \hat{\rho} - \hat{\sigma} \rVert_1$

Note that ρ and σ must be either Ket or Operator.

fidelity(ρ::QuantumObject, σ::QuantumObject)

Calculate the fidelity of two QuantumObject: $F(\hat{\rho}, \hat{\sigma}) = \textrm{Tr} \sqrt{\sqrt{\hat{\rho}} \hat{\sigma} \sqrt{\hat{\rho}}}$

Here, the definition is from Nielsen & Chuang, "Quantum Computation and Quantum Information". It is the square root of the fidelity defined in R. Jozsa, Journal of Modern Optics, 41:12, 2315 (1994).

Note that ρ and σ must be either Ket or Operator.

wigner(state::QuantumObject, xvec::AbstractVector, yvec::AbstractVector; g::Real=√2,
    solver::WignerSolver=WignerLaguerre())

Generates the Wigner quasipropability distribution of state at points xvec + 1im * yvec. The g parameter is a scaling factor related to the value of $\hbar$ in the commutation relation $[x, y] = i \hbar$ via $\hbar=2/g^2$ giving the default value $\hbar=1$.

The solver parameter can be either WignerLaguerre() or WignerClenshaw(). The former uses the Laguerre polynomial expansion of the Wigner function, while the latter uses the Clenshaw algorithm. The Laguerre expansion is faster for sparse matrices, while the Clenshaw algorithm is faster for dense matrices. The WignerLaguerre solver has an optional parallel parameter which defaults to true and uses multithreading to speed up the calculation.

negativity(ρ::QuantumObject, subsys::Int; logarithmic::Bool=false)

Compute the negativity $N(\hat{\rho}) = \frac{\Vert \hat{\rho}^{\Gamma}\Vert_1 - 1}{2}$ where $\hat{\rho}^{\Gamma}$ is the partial transpose of $\hat{\rho}$ with respect to the subsystem, and $\Vert \hat{X} \Vert_1=\textrm{Tr}\sqrt{\hat{X}^\dagger \hat{X}}$ is the trace norm.

Arguments

• ρ::QuantumObject: The density matrix (ρ.type must be OperatorQuantumObject).
• subsys::Int: an index that indicates which subsystem to compute the negativity for.
• logarithmic::Bool: choose whether to calculate logarithmic negativity or not. Default as false

Returns

• N::Real: The value of negativity.

Examples

julia> Ψ = bell_state(0, 0)
Quantum Object:   type=Ket   dims=[2, 2]   size=(4,)
4-element Vector{ComplexF64}:
 0.7071067811865475 + 0.0im
                0.0 + 0.0im
                0.0 + 0.0im
 0.7071067811865475 + 0.0im

julia> ρ = ket2dm(Ψ)
Quantum Object:   type=Operator   dims=[2, 2]   size=(4, 4)   ishermitian=true
4×4 Matrix{ComplexF64}:
 0.5+0.0im  0.0+0.0im  0.0+0.0im  0.5+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  0.0+0.0im
 0.5+0.0im  0.0+0.0im  0.0+0.0im  0.5+0.0im

julia> negativity(ρ, 2)
0.4999999999999998

AbstractLinearMap{T, TS}

Represents a general linear map with element type T and size TS.

Overview

A linear map is a transformation L that satisfies:

• Additivity: math L(u + v) = L(u) + L(v)
• Homogeneity: math L(cu) = cL(u)

It is typically represented as a matrix with dimensions given by size, and this abtract type helps to define this map when the matrix is not explicitly available.

Methods

• Base.eltype(A): Returns the element type T.
• Base.size(A): Returns the size A.size.
• Base.size(A, i): Returns the i-th dimension.

Example

As an example, we now define the linear map used in the eigsolve_al function for Arnoldi-Lindblad eigenvalue solver:

struct ArnoldiLindbladIntegratorMap{T,TS,TI} <: AbstractLinearMap{T,TS}
    elty::Type{T}
    size::TS
    integrator::TI
end

function LinearAlgebra.mul!(y::AbstractVector, A::ArnoldiLindbladIntegratorMap, x::AbstractVector)
    reinit!(A.integrator, x)
    solve!(A.integrator)
    return copyto!(y, A.integrator.u)
end

where integrator is the ODE integrator for the time-evolution. In this way, we can diagonalize this linear map using the eigsolve function.

QuantumToolbox.versioninfo(io::IO=stdout)

Command line output of information on QuantumToolbox, dependencies, and system information, same as QuantumToolbox.about.

QuantumToolbox.about(io::IO=stdout)

Command line output of information on QuantumToolbox, dependencies, and system information, same as QuantumToolbox.versioninfo.

gaussian(x::Number, μ::Number, σ::Number)

Returns the gaussian function $\exp \left[- \frac{(x - \mu)^2}{2 \sigma^2} \right]$, where $\mu$ and $\sigma^2$ are the mean and the variance respectively.

n_th(ω::Number, T::Real)

Gives the mean number of excitations in a mode with frequency ω at temperature T: $n_{\rm th} (\omega, T) = \frac{1}{e^{\omega/T} - 1}$

row_major_reshape(Q::AbstractArray, shapes...)

Reshapes Q in the row-major order, as numpy.

meshgrid(x::AbstractVector, y::AbstractVector)

Equivalent to numpy meshgrid.

_calculate_expectation!(p,z,B,idx) where T
Calculates the expectation values and function values of the operators and functions in p.e_ops and p.f_ops, respectively, and stores them in p.expvals and p.funvals.
The function is called by the callback _save_affect_lr_mesolve!.

Parameters
----------
p : NamedTuple
    The parameters of the problem.
z : AbstractMatrix{T}
    The z matrix.
B : AbstractMatrix{T}
    The B matrix.
idx : Integer
    The index of the current time step.

_adjM_condition_variational(u, t, integrator) where T
Condition for the dynamical rank adjustment based on the leakage out of the low-rank manifold.

Parameters
----------
u : AbstractVector{T}
    The current state of the system.
t : Real
    The current time.
integrator : ODEIntegrator
    The integrator of the problem.

_adjM_affect!(integrator)
Affect function for the dynamical rank adjustment. It increases the rank of the low-rank manifold by one, and updates the matrices accordingly.
If Δt>0, it rewinds the integrator to the previous time step.

Parameters
----------
integrator : ODEIntegrator
    The integrator of the problem.

_adjM_condition_ratio(u, t, integrator) where T
Condition for the dynamical rank adjustment based on the ratio between the smallest and largest eigenvalues of the density matrix.
The spectrum of the density matrix is calculated efficiently using the properties of the SVD decomposition of the matrix.

Parameters
----------
u : AbstractVector{T}
    The current state of the system.
t : Real
    The current time.
integrator : ODEIntegrator
    The integrator of the problem.

_pinv!(A, T1, T2; atol::Real=0.0, rtol::Real=(eps(real(float(oneunit(T))))*min(size(A)...))*iszero(atol)) where T
Computes the pseudo-inverse of a matrix A, and stores it in T1. If T2 is provided, it is used as a temporary matrix. 
The algorithm is based on the SVD decomposition of A, and is taken from the Julia package LinearAlgebra.
The difference with respect to the original function is that the cutoff is done with a smooth function instead of a step function.

Parameters
----------
A : AbstractMatrix{T}
    The matrix to be inverted.
T1 : AbstractMatrix{T}
T2 : AbstractMatrix{T}
    Temporary matrices used in the calculation.
atol : Real
    Absolute tolerance for the calculation of the pseudo-inverse.   
rtol : Real
    Relative tolerance for the calculation of the pseudo-inverse.

dBdz!(du, u, p, t) where T
Dynamical evolution equations for the low-rank manifold. The function is called by the ODEProblem.

Parameters
----------
du : AbstractVector{T}
    The derivative of the state of the system.
u : AbstractVector{T}
    The current state of the system.
p : NamedTuple
    The parameters of the problem.
t : Real
    The current time.