API

abstract type AbstractQuantumObject{DataType,ObjType,N}

Abstract type for all quantum objects like QuantumObject and QuantumObjectEvolution.

Example

julia> sigmax() isa AbstractQuantumObject
true

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 QuantumObjectEvolution{DT<:AbstractSciMLOperator,ObjType<:QuantumObjectType,N} <: AbstractQuantumObject
    data::DT
    type::ObjType
    dims::SVector{N,Int}
end

Julia struct representing any time-dependent quantum object. The data field is a AbstractSciMLOperator object that represents the time-dependent quantum object. It can be seen as

\[\hat{O}(t) = \sum_{i} c_i(p, t) \hat{O}_i\]

where $c_i(p, t)$ is a function that depends on the parameters p and time t, and $\hat{O}_i$ are the operators that form the quantum object. The type field is the type of the quantum object, and the dims field is the dimensions of the quantum object. For more information about type and dims, see QuantumObject. For more information about AbstractSciMLOperator, see the SciML documentation.

Examples

This operator can be initialized in the same way as the QuTiP QobjEvo object. For example

julia> a = tensor(destroy(10), qeye(2))
Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=false
20×20 SparseMatrixCSC{ComplexF64, Int64} with 18 stored entries:
⎡⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⎦

julia> σm = tensor(qeye(10), sigmam())
Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=false
20×20 SparseMatrixCSC{ComplexF64, Int64} with 10 stored entries:
⎡⠂⡀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠀⠂⡀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠂⡀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠂⡀⠀⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠂⡀⎦

julia> coef1(p, t) = exp(-1im * t)
coef1 (generic function with 1 method)

julia> coef2(p, t) = sin(t)
coef2 (generic function with 1 method)

julia> op1 = QuantumObjectEvolution(((a, coef1), (σm, coef2)))
Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=true
(ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20) + ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20))

We can also concretize the operator at a specific time t

julia> op1(0.1)
Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=false
20×20 SparseMatrixCSC{ComplexF64, Int64} with 28 stored entries:
⎡⠂⡑⢄⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠀⠂⡑⢄⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠂⡑⢄⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠂⡑⢄⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠂⡑⎦

It also supports parameter-dependent time evolution

julia> coef1(p, t) = exp(-1im * p.ω1 * t)
coef1 (generic function with 1 method)

julia> coef2(p, t) = sin(p.ω2 * t)
coef2 (generic function with 1 method)

julia> op1 = QuantumObjectEvolution(((a, coef1), (σm, coef2)))
Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=true
(ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20) + ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20))

julia> p = (ω1 = 1.0, ω2 = 0.5)
(ω1 = 1.0, ω2 = 0.5)

julia> op1(p, 0.1)
Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=false
20×20 SparseMatrixCSC{ComplexF64, Int64} with 28 stored entries:
⎡⠂⡑⢄⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠀⠂⡑⢄⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠂⡑⢄⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠂⡑⢄⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠂⡑⎦

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

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

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

eltype(A::AbstractQuantumObject)

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

length(A::AbstractQuantumObject)

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

isbra(A)

Checks if the QuantumObject A is a BraQuantumObject. Default case returns false for any other inputs.

isket(A)

Checks if the QuantumObject A is a KetQuantumObject. Default case returns false for any other inputs.

isoper(A)

Checks if the AbstractQuantumObject A is a OperatorQuantumObject. Default case returns false for any other inputs.

isoperbra(A)

Checks if the QuantumObject A is a OperatorBraQuantumObject. Default case returns false for any other inputs.

isoperket(A)

Checks if the QuantumObject A is a OperatorKetQuantumObject. Default case returns false for any other inputs.

issuper(A)

Checks if the AbstractQuantumObject A is a SuperOperatorQuantumObject. Default case returns false for any other inputs.

ishermitian(A::AbstractQuantumObject)

Test whether the AbstractQuantumObject is Hermitian.

issymmetric(A::AbstractQuantumObject)

Test whether the AbstractQuantumObject is symmetric.

isposdef(A::AbstractQuantumObject)

Test whether the AbstractQuantumObject 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.

isconstant(A::AbstractQuantumObject)

Test whether the AbstractQuantumObject A is constant in time. For a QuantumObject, this function returns true, while for a QuantumObjectEvolution, this function returns true if the operator is contant in time.

conj(A::AbstractQuantumObject)

Return the element-wise complex conjugation of the AbstractQuantumObject.

transpose(A::AbstractQuantumObject)

Lazy matrix transpose of the AbstractQuantumObject.

A'
adjoint(A::AbstractQuantumObject)

Lazy adjoint (conjugate transposition) of the AbstractQuantumObject

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::AbstractQuantumObject j::QuantumObject)

Compute the generalized dot product dot(i, A*j) between a AbstractQuantumObject and two QuantumObject (i and j), namely $\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::AbstractQuantumObject)

Matrix inverse of the AbstractQuantumObject. If A is a QuantumObjectEvolution, the inverse is computed at the last computed time.

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::AbstractQuantumObject)

Returns the data of a AbstractQuantumObject.

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::Union{AbstractQuantumObject{DT1,HOpType},Tuple},
    T::Real,
    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
    alg::OrdinaryDiffEqAlgorithm = Tsit5(),
    params::NamedTuple = NamedTuple(),
    ρ0::AbstractMatrix = rand_dm(prod(H.dims)).data,
    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. It can be a QuantumObject, a QuantumObjectEvolution, or a tuple of the form supported by mesolve.
• 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. Default is Tsit5().
• params: A NamedTuple containing the parameters of the system.
• ρ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::AbstractQuantumObject, ψ::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::AbstractQuantumObject, B::AbstractQuantumObject, ...)

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::AbstractQuantumObject, 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::AbstractQuantumObject, 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, spost, and sprepost.

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

Generate QuantumObject

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

QobjEvo(op_func_list::Union{Tuple,AbstractQuantumObject}, α::Union{Nothing,Number}=nothing; type::Union{Nothing, QuantumObjectType}=nothing, f::Function=identity)

Generate QuantumObjectEvolution.

Note that this functions is same as QuantumObjectEvolution(op_func_list). If α is provided, all the operators in op_func_list will be pre-multiplied by α. The type parameter is used to specify the type of the QuantumObject, either Operator or SuperOperator. The f parameter is used to pre-apply a function to the operators before converting them to SciML operators.

Arguments

• op_func_list::Union{Tuple,AbstractQuantumObject}: A tuple of tuples or operators.
• α::Union{Nothing,Number}=nothing: A scalar to pre-multiply the operators.
• f::Function=identity: A function to pre-apply to the operators.

Examples

This operator can be initialized in the same way as the QuTiP QobjEvo object. For example

julia> a = tensor(destroy(10), qeye(2))
Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=false
20×20 SparseMatrixCSC{ComplexF64, Int64} with 18 stored entries:
⎡⠀⠑⢄⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠑⢄⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠀⠑⎦

julia> σm = tensor(qeye(10), sigmam())
Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=false
20×20 SparseMatrixCSC{ComplexF64, Int64} with 10 stored entries:
⎡⠂⡀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠀⠂⡀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠂⡀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠂⡀⠀⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠂⡀⎦

julia> coef1(p, t) = exp(-1im * t)
coef1 (generic function with 1 method)

julia> coef2(p, t) = sin(t)
coef2 (generic function with 1 method)

julia> op1 = QobjEvo(((a, coef1), (σm, coef2)))
Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=true
(ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20) + ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20))

We can also concretize the operator at a specific time t

julia> op1(0.1)
Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=false
20×20 SparseMatrixCSC{ComplexF64, Int64} with 28 stored entries:
⎡⠂⡑⢄⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠀⠂⡑⢄⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠂⡑⢄⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠂⡑⢄⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠂⡑⎦

It also supports parameter-dependent time evolution

julia> coef1(p, t) = exp(-1im * p.ω1 * t)
coef1 (generic function with 1 method)

julia> coef2(p, t) = sin(p.ω2 * t)
coef2 (generic function with 1 method)

julia> op1 = QobjEvo(((a, coef1), (σm, coef2)))
Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=true
(ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20) + ScalarOperator(0.0 + 0.0im) * MatrixOperator(20 × 20))

julia> p = (ω1 = 1.0, ω2 = 0.5)
(ω1 = 1.0, ω2 = 0.5)

julia> op1(p, 0.1)
Quantum Object:   type=Operator   dims=[10, 2]   size=(20, 20)   ishermitian=false
20×20 SparseMatrixCSC{ComplexF64, Int64} with 28 stored entries:
⎡⠂⡑⢄⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠀⠀⠂⡑⢄⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠂⡑⢄⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠂⡑⢄⠀⎥
⎣⠀⠀⠀⠀⠀⠀⠀⠀⠂⡑⎦

QobjEvo(data::AbstractSciMLOperator; type::QuantumObjectType = Operator, dims = nothing)

Synonym of QuantumObjectEvolution object from a SciMLOperator. See the documentation for QuantumObjectEvolution for more information.

shape(A::AbstractQuantumObject)

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

Note that this function is same as size(A).

isherm(A::AbstractQuantumObject)

Test whether the AbstractQuantumObject is Hermitian.

Note that this functions is same as ishermitian(A).

trans(A::AbstractQuantumObject)

Lazy matrix transpose of the AbstractQuantumObject.

Note that this function is same as transpose(A).

dag(A::AbstractQuantumObject)

Lazy adjoint (conjugate transposition) of the AbstractQuantumObject

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::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
    ψ0::QuantumObject{DT2,KetQuantumObject},
    tlist::AbstractVector;
    e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
    params::NamedTuple = NamedTuple(),
    progress_bar::Union{Val,Bool} = Val(true),
    kwargs...,
)

Generate 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: Hamiltonian of the system $\hat{H}$. It can be either a QuantumObject, a QuantumObjectEvolution, or a Tuple of operator-function pairs.
• ψ0: Initial state of the system $|\psi(0)\rangle$.
• tlist: List of times at which to save either the state or the expectation values of the system.
• e_ops: List of operators for which to calculate expectation values. It can be either a Vector or a Tuple.
• params: NamedTuple of parameters to pass to the solver.
• progress_bar: 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 empty, the default value of saveat=tlist (saving the states corresponding to tlist), otherwise, saveat=[tlist[end]] (only save the final state). 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 kwargs please refer to DifferentialEquations.jl (Keyword Arguments)

Returns

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

mesolveProblem(
    H::Union{AbstractQuantumObject{DT1,HOpType},Tuple},
    ψ0::QuantumObject{DT2,StateOpType},
    tlist,
    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
    e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
    params::NamedTuple = NamedTuple(),
    progress_bar::Union{Val,Bool} = Val(true),
    kwargs...,
)

Generate 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: Hamiltonian of the system $\hat{H}$. It can be either a QuantumObject, a QuantumObjectEvolution, or a Tuple of operator-function pairs.
• ψ0: Initial state of the system $|\psi(0)\rangle$.
• tlist: List of times at which to save either the state or the expectation values of the system.
• c_ops: List of collapse operators $\{\hat{C}_n\}_n$. It can be either a Vector or a Tuple.
• e_ops: List of operators for which to calculate expectation values. It can be either a Vector or a Tuple.
• params: NamedTuple of parameters to pass to the solver.
• progress_bar: 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 empty, the default value of saveat=tlist (saving the states corresponding to tlist), otherwise, saveat=[tlist[end]] (only save the final state). 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 kwargs please refer to DifferentialEquations.jl (Keyword Arguments)

Returns

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

mcsolveProblem(
    H::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
    ψ0::QuantumObject{DT2,KetQuantumObject},
    tlist::AbstractVector,
    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
    e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
    params::NamedTuple = NamedTuple(),
    rng::AbstractRNG = default_rng(),
    jump_callback::TJC = ContinuousLindbladJumpCallback(),
    kwargs...,
)

Generate 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: Hamiltonian of the system $\hat{H}$. It can be either a QuantumObject, a QuantumObjectEvolution, or a Tuple of operator-function pairs.
• ψ0: Initial state of the system $|\psi(0)\rangle$.
• tlist: List of times at which to save either the state or the expectation values of the system.
• c_ops: List of collapse operators $\{\hat{C}_n\}_n$. It can be either a Vector or a Tuple.
• e_ops: List of operators for which to calculate expectation values. It can be either a Vector or a Tuple.
• params: NamedTuple of parameters to pass to the solver.
• rng: Random number generator for reproducibility.
• jump_callback: The Jump Callback type: Discrete or Continuous. The default is ContinuousLindbladJumpCallback(), which is more precise.
• 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 empty, the default value of saveat=tlist (saving the states corresponding to tlist), otherwise, saveat=[tlist[end]] (only save the final state). 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 kwargs please refer to DifferentialEquations.jl (Keyword Arguments)

Returns

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

mcsolveEnsembleProblem(
    H::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
    ψ0::QuantumObject{DT2,KetQuantumObject},
    tlist::AbstractVector,
    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
    e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
    params::NamedTuple = NamedTuple(),
    rng::AbstractRNG = default_rng(),
    ntraj::Int = 1,
    ensemble_method = EnsembleThreads(),
    jump_callback::TJC = ContinuousLindbladJumpCallback(),
    prob_func::Function = _mcsolve_prob_func,
    output_func::Function = _mcsolve_dispatch_output_func(ensemble_method),
    progress_bar::Union{Val,Bool} = Val(true),
    kwargs...,
)

Generate 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: Hamiltonian of the system $\hat{H}$. It can be either a QuantumObject, a QuantumObjectEvolution, or a Tuple of operator-function pairs.
• ψ0: Initial state of the system $|\psi(0)\rangle$.
• tlist: List of times at which to save either the state or the expectation values of the system.
• c_ops: List of collapse operators $\{\hat{C}_n\}_n$. It can be either a Vector or a Tuple.
• e_ops: List of operators for which to calculate expectation values. It can be either a Vector or a Tuple.
• params: NamedTuple of parameters to pass to the solver.
• rng: Random number generator for reproducibility.
• ntraj: Number of trajectories to use.
• ensemble_method: Ensemble method to use. Default to EnsembleThreads().
• jump_callback: The Jump Callback type: Discrete or Continuous. The default is ContinuousLindbladJumpCallback(), which is more precise.
• prob_func: Function to use for generating the ODEProblem.
• output_func: Function to use for generating the output of a single trajectory.
• progress_bar: 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 empty, the default value of saveat=tlist (saving the states corresponding to tlist), otherwise, saveat=[tlist[end]] (only save the final state). 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 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::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
    ψ0::QuantumObject{DT2,KetQuantumObject},
    tlist::AbstractVector,
    sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
    e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
    params::NamedTuple = NamedTuple(),
    rng::AbstractRNG = default_rng(),
    progress_bar::Union{Val,Bool} = Val(true),
    kwargs...,
)

Generate 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

\[K = \hat{H} + i \sum_n \left(\frac{e_j} C_n - \frac{1}{2} \sum_{j} C_n^\dagger C_n - \frac{e_j^2}{8}\right),\]

\[M_n = C_n - \frac{e_n}{2},\]

and

\[e_n = \langle C_n + 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: Hamiltonian of the system $\hat{H}$. It can be either a QuantumObject, a QuantumObjectEvolution, or a Tuple of operator-function pairs.
• ψ0: Initial state of the system $|\psi(0)\rangle$.
• tlist: List of times at which to save either the state or the expectation values of the system.
• sc_ops: List of collapse operators $\{\hat{C}_n\}_n$. It can be either a Vector or a Tuple.
• e_ops: List of operators for which to calculate expectation values. It can be either a Vector or a Tuple.
• params: NamedTuple of parameters to pass to the solver.
• rng: Random number generator for reproducibility.
• progress_bar: 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 empty, the default value of saveat=tlist (saving the states corresponding to tlist), otherwise, saveat=[tlist[end]] (only save the final state). 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 kwargs please refer to DifferentialEquations.jl (Keyword Arguments)

Returns

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

ssesolveEnsembleProblem(
    H::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
    ψ0::QuantumObject{DT2,KetQuantumObject},
    tlist::AbstractVector,
    sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
    e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
    params::NamedTuple = NamedTuple(),
    rng::AbstractRNG = default_rng(),
    ntraj::Int = 1,
    ensemble_method = EnsembleThreads(),
    prob_func::Function = _ssesolve_prob_func,
    output_func::Function = _ssesolve_dispatch_output_func(ensemble_method),
    progress_bar::Union{Val,Bool} = Val(true),
    kwargs...,
)

Generate 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

\[K = \hat{H} + i \sum_n \left(\frac{e_j} C_n - \frac{1}{2} \sum_{j} C_n^\dagger C_n - \frac{e_j^2}{8}\right),\]

\[M_n = C_n - \frac{e_n}{2},\]

and

\[e_n = \langle C_n + 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: Hamiltonian of the system $\hat{H}$. It can be either a QuantumObject, a QuantumObjectEvolution, or a Tuple of operator-function pairs.
• ψ0: Initial state of the system $|\psi(0)\rangle$.
• tlist: List of times at which to save either the state or the expectation values of the system.
• sc_ops: List of collapse operators $\{\hat{C}_n\}_n$. It can be either a Vector or a Tuple.
• e_ops: List of operators for which to calculate expectation values. It can be either a Vector or a Tuple.
• params: NamedTuple of parameters to pass to the solver.
• rng: Random number generator for reproducibility.
• ntraj: Number of trajectories to use.
• ensemble_method: Ensemble method to use. Default to EnsembleThreads().
• jump_callback: The Jump Callback type: Discrete or Continuous. The default is ContinuousLindbladJumpCallback(), which is more precise.
• prob_func: Function to use for generating the ODEProblem.
• output_func: Function to use for generating the output of a single trajectory.
• progress_bar: 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 empty, the default value of saveat=tlist (saving the states corresponding to tlist), otherwise, saveat=[tlist[end]] (only save the final state). 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 kwargs please refer to DifferentialEquations.jl (Keyword Arguments)

Returns

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

sesolve(
    H::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
    ψ0::QuantumObject{DT2,KetQuantumObject},
    tlist::AbstractVector;
    alg::OrdinaryDiffEqAlgorithm = Tsit5(),
    e_ops::Union{Nothing,AbstractVector,Tuple} = 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: Hamiltonian of the system $\hat{H}$. It can be either a QuantumObject, a QuantumObjectEvolution, or a Tuple of operator-function pairs.
• ψ0: Initial state of the system $|\psi(0)\rangle$.
• tlist: List of times at which to save either the state or the expectation values of the system.
• alg: The algorithm for the ODE solver. The default is Tsit5().
• e_ops: List of operators for which to calculate expectation values. It can be either a Vector or a Tuple.
• params: NamedTuple of parameters to pass to the solver.
• progress_bar: 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 empty, the default value of saveat=tlist (saving the states corresponding to tlist), otherwise, saveat=[tlist[end]] (only save the final state). 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::Union{AbstractQuantumObject{DT1,HOpType},Tuple},
    ψ0::QuantumObject{DT2,StateOpType},
    tlist::AbstractVector,
    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
    alg::OrdinaryDiffEqAlgorithm = Tsit5(),
    e_ops::Union{Nothing,AbstractVector,Tuple} = 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: Hamiltonian of the system $\hat{H}$. It can be either a QuantumObject, a QuantumObjectEvolution, or a Tuple of operator-function pairs.
• ψ0: Initial state of the system $|\psi(0)\rangle$.
• tlist: List of times at which to save either the state or the expectation values of the system.
• c_ops: List of collapse operators $\{\hat{C}_n\}_n$. It can be either a Vector or a Tuple.
• alg: The algorithm for the ODE solver. The default value is Tsit5().
• e_ops: List of operators for which to calculate expectation values. It can be either a Vector or a Tuple.
• params: NamedTuple of parameters to pass to the solver.
• progress_bar: 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 empty, the default value of saveat=tlist (saving the states corresponding to tlist), otherwise, saveat=[tlist[end]] (only save the final state). 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::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
    ψ0::QuantumObject{DT2,KetQuantumObject},
    tlist::AbstractVector,
    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
    alg::OrdinaryDiffEqAlgorithm = Tsit5(),
    e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
    params::NamedTuple = NamedTuple(),
    rng::AbstractRNG = default_rng(),
    ntraj::Int = 1,
    ensemble_method = EnsembleThreads(),
    jump_callback::TJC = ContinuousLindbladJumpCallback(),
    prob_func::Function = _mcsolve_prob_func,
    output_func::Function = _mcsolve_dispatch_output_func(ensemble_method),
    progress_bar::Union{Val,Bool} = Val(true),
    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: Hamiltonian of the system $\hat{H}$. It can be either a QuantumObject, a QuantumObjectEvolution, or a Tuple of operator-function pairs.
• ψ0: Initial state of the system $|\psi(0)\rangle$.
• tlist: List of times at which to save either the state or the expectation values of the system.
• c_ops: List of collapse operators $\{\hat{C}_n\}_n$. It can be either a Vector or a Tuple.
• alg: The algorithm to use for the ODE solver. Default to Tsit5().
• e_ops: List of operators for which to calculate expectation values. It can be either a Vector or a Tuple.
• params: NamedTuple of parameters to pass to the solver.
• rng: Random number generator for reproducibility.
• ntraj: Number of trajectories to use.
• ensemble_method: Ensemble method to use. Default to EnsembleThreads().
• jump_callback: The Jump Callback type: Discrete or Continuous. The default is ContinuousLindbladJumpCallback(), which is more precise.
• prob_func: Function to use for generating the ODEProblem.
• output_func: Function to use for generating the output of a single trajectory.
• progress_bar: Whether to show the progress bar. Using non-Val types might lead to type instabilities.
• kwargs: The keyword arguments for the ODEProblem.

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 empty, the default value of saveat=tlist (saving the states corresponding to tlist), otherwise, saveat=[tlist[end]] (only save the final state). 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::Union{AbstractQuantumObject{DT1,OperatorQuantumObject},Tuple},
    ψ0::QuantumObject{DT2,KetQuantumObject},
    tlist::AbstractVector,
    sc_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
    alg::StochasticDiffEqAlgorithm = SRA1(),
    e_ops::Union{Nothing,AbstractVector,Tuple} = nothing,
    params::NamedTuple = NamedTuple(),
    rng::AbstractRNG = default_rng(),
    ntraj::Int = 1,
    ensemble_method = EnsembleThreads(),
    prob_func::Function = _ssesolve_prob_func,
    output_func::Function = _ssesolve_dispatch_output_func(ensemble_method),
    progress_bar::Union{Val,Bool} = Val(true),
    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

\[K = \hat{H} + i \sum_n \left(\frac{e_j} C_n - \frac{1}{2} \sum_{j} C_n^\dagger C_n - \frac{e_j^2}{8}\right),\]

\[M_n = C_n - \frac{e_n}{2},\]

and

\[e_n = \langle C_n + 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: Hamiltonian of the system $\hat{H}$. It can be either a QuantumObject, a QuantumObjectEvolution, or a Tuple of operator-function pairs.
• ψ0: Initial state of the system $|\psi(0)\rangle$.
• tlist: List of times at which to save either the state or the expectation values of the system.
• sc_ops: List of collapse operators $\{\hat{C}_n\}_n$. It can be either a Vector or a Tuple.
• alg: The algorithm to use for the stochastic differential equation. Default is SRA1().
• e_ops: List of operators for which to calculate expectation values. It can be either a Vector or a Tuple.
• params: NamedTuple of parameters to pass to the solver.
• rng: Random number generator for reproducibility.
• ntraj: Number of trajectories to use.
• ensemble_method: Ensemble method to use. Default to EnsembleThreads().
• prob_func: Function to use for generating the ODEProblem.
• output_func: Function to use for generating the output of a single trajectory.
• progress_bar: 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 empty, the default value of saveat=tlist (saving the states corresponding to tlist), otherwise, saveat=[tlist[end]] (only save the final state). 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,
    tlist::AbstractVector, c_ops::Function, maxdims::AbstractVector,
    dfd_params::NamedTuple=NamedTuple();
    alg::OrdinaryDiffEqAlgorithm=Tsit5(),
    e_ops::Function=(dim_list) -> Vector{Vector{T1}}([]),
    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)

liouvillian(H::AbstractQuantumObject, c_ops::Union{Nothing,AbstractVector,Tuple}=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{DT,OpType},
    c_ops::Union{Nothing,AbstractVector,Tuple} = nothing;
    solver::SteadyStateSolver = SteadyStateDirectSolver(),
    kwargs...,
)

Solve the stationary state based on different solvers.

Parameters

• H: The Hamiltonian or the Liouvillian of the system.
• c_ops: The list of the collapse operators.
• solver: see documentation Solving for Steady-State Solutions for different solvers.
• kwargs: The keyword arguments for the solver.

steadystate(
    H::QuantumObject{DT1,HOpType},
    ψ0::QuantumObject{DT2,StateOpType},
    tmax::Real = Inf,
    c_ops::Union{Nothing,AbstractVector,Tuple} = 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: The Hamiltonian or the Liouvillian of the system.
• ψ0: The initial state of the system.
• tmax=Inf: The final time step for the steady state problem.
• c_ops=nothing: The list of the collapse operators.
• solver: see SteadyStateODESolver for more details.
• reltol=1.0e-8: Relative tolerance in steady state terminate condition and solver adaptive timestepping.
• abstol=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,Tuple} = 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::Union{Nothing,AbstractVector} = nothing: 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

dsf_mesolve(H::Function,
    ψ0::QuantumObject,
    tlist::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}([]),
    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,
    tlist::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}([]),
    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))),
    progress_bar::Union{Bool,Val} = Val(true)
    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)

struct TimeEvolutionLRSol

A structure storing the results and some information from solving low-rank master equation 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.
• fexpect::Matrix: The function values at each time point.
• 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.
• z::Vector{QuantumObject}: The z` matrix of the low-rank algorithm at each time point.
• B::Vector{QuantumObject}: The B matrix of the low-rank algorithm at each time point.

lr_mesolveProblem(
    H::QuantumObject{DT1,OperatorQuantumObject},
    z::AbstractArray{T2,2},
    B::AbstractArray{T2,2},
    tlist::AbstractVector,
    c_ops::Union{AbstractVector,Tuple}=();
    e_ops::Union{AbstractVector,Tuple}=(),
    f_ops::Union{AbstractVector,Tuple}=(),
    opt::NamedTuple = lr_mesolve_options_default,
    kwargs...,
)

Formulates the ODEproblem for the low-rank time evolution of the system. The function is called by lr_mesolve. For more information about the low-rank master equation, see (Gravina and Savona, 2024).

Arguments

• H::QuantumObject: The Hamiltonian of the system.
• z::AbstractArray: The initial z matrix of the low-rank algorithm.
• B::AbstractArray: The initial B matrix of the low-rank algorithm.
• tlist::AbstractVector: The time steps at which the expectation values and function values are calculated.
• c_ops::Union{AbstractVector,Tuple}: The list of the collapse operators.
• e_ops::Union{AbstractVector,Tuple}: The list of the operators for which the expectation values are calculated.
• f_ops::Union{AbstractVector,Tuple}: The list of the functions for which the function values are calculated.
• opt::NamedTuple: The options of the low-rank master equation.
• kwargs: Additional keyword arguments.

lr_mesolve(
    H::QuantumObject{DT1,OperatorQuantumObject},
    z::AbstractArray{T2,2},
    B::AbstractArray{T2,2},
    tlist::AbstractVector,
    c_ops::Union{AbstractVector,Tuple}=();
    e_ops::Union{AbstractVector,Tuple}=(),
    f_ops::Union{AbstractVector,Tuple}=(),
    opt::NamedTuple = lr_mesolve_options_default,
    kwargs...,
)

Time evolution of an open quantum system using the low-rank master equation. For more information about the low-rank master equation, see (Gravina and Savona, 2024).

Arguments

• H::QuantumObject: The Hamiltonian of the system.
• z::AbstractArray: The initial z matrix of the low-rank algorithm.
• B::AbstractArray: The initial B matrix of the low-rank algorithm.
• tlist::AbstractVector: The time steps at which the expectation values and function values are calculated.
• c_ops::Union{AbstractVector,Tuple}: The list of the collapse operators.
• e_ops::Union{AbstractVector,Tuple}: The list of the operators for which the expectation values are calculated.
• f_ops::Union{AbstractVector,Tuple}: The list of the functions for which the function values are calculated.
• opt::NamedTuple: The options of the low-rank master equation.
• kwargs: Additional keyword arguments.

correlation_3op_2t(H::QuantumObject,
    ψ0::QuantumObject,
    t_l::AbstractVector,
    τ_l::AbstractVector,
    A::QuantumObject,
    B::QuantumObject,
    C::QuantumObject,
    c_ops::Union{Nothing,AbstractVector,Tuple}=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,Tuple}=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,Tuple}=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,Tuple}=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.

Lattice

A Julia constructor for a lattice object. The lattice object is used to define the geometry of the lattice. Nx and Ny are the number of sites in the x and y directions, respectively. N is the total number of sites. lin_idx is a LinearIndices object and car_idx is a CartesianIndices object, and they are used to efficiently select sites on the lattice.

SingleSiteOperator(O::QuantumObject, i::Integer, N::Integer)

A Julia constructor for a single-site operator. s is the operator acting on the site. i is the site index, and N is the total number of sites. The function returns a QuantumObject given by $\\mathbb{1}^{\\otimes (i - 1)} \\otimes \hat{O} \\otimes \\mathbb{1}^{\\otimes (N - i)}$.

DissipativeIsing(Jx::Real, Jy::Real, Jz::Real, hx::Real, hy::Real, hz::Real, γ::Real, latt::Lattice; boundary_condition::Union{Symbol, Val} = Val(:periodic_bc), order::Integer = 1)

A Julia constructor for a dissipative Ising model. The function returns the Hamiltonian

\[\hat{H} = \frac{J_x}{2} \sum_{\langle i, j \rangle} \hat{\sigma}_i^x \hat{\sigma}_j^x + \frac{J_y}{2} \sum_{\langle i, j \rangle} \hat{\sigma}_i^y \hat{\sigma}_j^y + \frac{J_z}{2} \sum_{\langle i, j \rangle} \hat{\sigma}_i^z \hat{\sigma}_j^z + h_x \sum_i \hat{\sigma}_i^x\]

and the collapse operators

\[\hat{c}_i = \sqrt{\gamma} \hat{\sigma}_i^-\]

Arguments

• Jx::Real: The coupling constant in the x-direction.
• Jy::Real: The coupling constant in the y-direction.
• Jz::Real: The coupling constant in the z-direction.
• hx::Real: The magnetic field in the x-direction.
• hy::Real: The magnetic field in the y-direction.
• hz::Real: The magnetic field in the z-direction.
• γ::Real: The local dissipation rate.
• latt::Lattice: A Lattice object that defines the geometry of the lattice.
• boundary_condition::Union{Symbol, Val}: The boundary conditions of the lattice. The possible inputs are periodic_bc and open_bc, for periodic or open boundary conditions, respectively. The default value is Val(:periodic_bc).
• order::Integer: The order of the nearest-neighbour sites. The default value is 1.

wigner(
    state::QuantumObject{DT,OpType},
    xvec::AbstractVector,
    yvec::AbstractVector;
    g::Real = √2,
    method::WignerSolver = WignerClenshaw(),
)

Generates the Wigner quasipropability distribution of state at points xvec + 1im * yvec in phase space. 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 method 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 method has an optional parallel parameter which defaults to true and uses multithreading to speed up the calculation.

Arguments

• state::QuantumObject: The quantum state for which the Wigner function is calculated. It can be either a KetQuantumObject, BraQuantumObject, or OperatorQuantumObject.
• xvec::AbstractVector: The x-coordinates of the phase space grid.
• yvec::AbstractVector: The y-coordinates of the phase space grid.
• g::Real: The scaling factor related to the value of $\hbar$ in the commutation relation $[x, y] = i \hbar$ via $\hbar=2/g^2$.
• method::WignerSolver: The method used to calculate the Wigner function. It can be either WignerLaguerre() or WignerClenshaw(), with WignerClenshaw() as default. The WignerLaguerre method has the optional parallel and tol parameters, with default values true and 1e-14, respectively.

Returns

• W::Matrix: The Wigner function of the state at the points xvec + 1im * yvec in phase space.

Example

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

julia> xvec = range(-5, 5, 200)
-5.0:0.05025125628140704:5.0

julia> wig = wigner(ψ, xvec, xvec)
200×200 Matrix{Float64}:
 2.63558e-21  4.30187e-21  6.98638e-21  1.12892e-20  1.81505e-20  …  1.50062e-20  9.28736e-21  5.71895e-21  3.50382e-21
 4.29467e-21  7.00905e-21  1.13816e-20  1.83891e-20  2.9562e-20      2.45173e-20  1.51752e-20  9.3454e-21   5.72614e-21
 6.96278e-21  1.13621e-20  1.8448e-20   2.98026e-20  4.79043e-20     3.98553e-20  2.46711e-20  1.51947e-20  9.31096e-21
 1.12314e-20  1.83256e-20  2.97505e-20  4.80558e-20  7.72344e-20     6.4463e-20   3.99074e-20  2.45808e-20  1.50639e-20
 1.80254e-20  2.94073e-20  4.77351e-20  7.70963e-20  1.23892e-19     1.0374e-19   6.42289e-20  3.95652e-20  2.42491e-20
 ⋮                                                                ⋱                                         
 1.80254e-20  2.94073e-20  4.77351e-20  7.70963e-20  1.23892e-19  …  1.0374e-19   6.42289e-20  3.95652e-20  2.42491e-20
 1.12314e-20  1.83256e-20  2.97505e-20  4.80558e-20  7.72344e-20     6.4463e-20   3.99074e-20  2.45808e-20  1.50639e-20
 6.96278e-21  1.13621e-20  1.8448e-20   2.98026e-20  4.79043e-20     3.98553e-20  2.46711e-20  1.51947e-20  9.31096e-21
 4.29467e-21  7.00905e-21  1.13816e-20  1.83891e-20  2.9562e-20      2.45173e-20  1.51752e-20  9.3454e-21   5.72614e-21
 2.63558e-21  4.30187e-21  6.98638e-21  1.12892e-20  1.81505e-20     1.50062e-20  9.28736e-21  5.71895e-21  3.50382e-21

or taking advantage of the parallel computation of the WignerLaguerre method

julia> wig = wigner(ρ, xvec, xvec, method=WignerLaguerre(parallel=true))
200×200 Matrix{Float64}:
 2.63558e-21  4.30187e-21  6.98638e-21  1.12892e-20  1.81505e-20  …  1.50062e-20  9.28736e-21  5.71895e-21  3.50382e-21
 4.29467e-21  7.00905e-21  1.13816e-20  1.83891e-20  2.9562e-20      2.45173e-20  1.51752e-20  9.3454e-21   5.72614e-21
 6.96278e-21  1.13621e-20  1.8448e-20   2.98026e-20  4.79043e-20     3.98553e-20  2.46711e-20  1.51947e-20  9.31096e-21
 1.12314e-20  1.83256e-20  2.97505e-20  4.80558e-20  7.72344e-20     6.4463e-20   3.99074e-20  2.45808e-20  1.50639e-20
 1.80254e-20  2.94073e-20  4.77351e-20  7.70963e-20  1.23892e-19     1.0374e-19   6.42289e-20  3.95652e-20  2.42491e-20
 ⋮                                                                ⋱                                         
 1.80254e-20  2.94073e-20  4.77351e-20  7.70963e-20  1.23892e-19  …  1.0374e-19   6.42289e-20  3.95652e-20  2.42491e-20
 1.12314e-20  1.83256e-20  2.97505e-20  4.80558e-20  7.72344e-20     6.4463e-20   3.99074e-20  2.45808e-20  1.50639e-20
 6.96278e-21  1.13621e-20  1.8448e-20   2.98026e-20  4.79043e-20     3.98553e-20  2.46711e-20  1.51947e-20  9.31096e-21
 4.29467e-21  7.00905e-21  1.13816e-20  1.83891e-20  2.9562e-20      2.45173e-20  1.51752e-20  9.3454e-21   5.72614e-21
 2.63558e-21  4.30187e-21  6.98638e-21  1.12892e-20  1.81505e-20     1.50062e-20  9.28736e-21  5.71895e-21  3.50382e-21

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_thermal(ω::Real, ω_th::Real)

Return the number of photons in thermal equilibrium for an harmonic oscillator mode with frequency $\omega$, at the temperature described by $\omega_{\textrm{th}} \equiv k_B T / \hbar$:

\[n(\omega, \omega_{\textrm{th}}) = \frac{1}{e^{\omega/\omega_{\textrm{th}}} - 1},\]

where $\hbar$ is the reduced Planck constant, and $k_B$ is the Boltzmann constant.

const PhysicalConstants

A NamedTuple which stores some constant values listed in CODATA recommended values of the fundamental physical constants: 2022

The current stored constants are:

• c : (exact) speed of light in vacuum with unit $[\textrm{m}\cdot\textrm{s}^{-1}]$
• G : Newtonian constant of gravitation with unit $[\textrm{m}^3\cdot\textrm{kg}^{−1}\cdot\textrm{s}^{−2}]$
• h : (exact) Planck constant with unit $[\textrm{J}\cdot\textrm{s}]$
• ħ : reduced Planck constant with unit $[\textrm{J}\cdot\textrm{s}]$
• e : (exact) elementary charge with unit $[\textrm{C}]$
• μ0 : vacuum magnetic permeability with unit $[\textrm{N}\cdot\textrm{A}^{-2}]$
• ϵ0 : vacuum electric permittivity with unit $[\textrm{F}\cdot\textrm{m}^{-1}]$
• k : (exact) Boltzmann constant with unit $[\textrm{J}\cdot\textrm{K}^{-1}]$
• NA : (exact) Avogadro constant with unit $[\textrm{mol}^{-1}]$

Examples

julia> PhysicalConstants.ħ
1.0545718176461565e-34

convert_unit(value::Real, unit1::Symbol, unit2::Symbol)

Convert the energy value from unit1 to unit2. The unit1 and unit2 can be either the following Symbol:

• :J : Joule
• :eV : electron volt
• :meV : milli-electron volt
• :MHz : Mega-Hertz multiplied by Planck constant $h$
• :GHz : Giga-Hertz multiplied by Planck constant $h$
• :K : Kelvin multiplied by Boltzmann constant $k$
• :mK : milli-Kelvin multiplied by Boltzmann constant $k$

Note that we use the values stored in PhysicalConstants to do the conversion.

Examples

julia> convert_unit(1, :eV, :J)
1.602176634e-19

julia> convert_unit(1, :GHz, :J)
6.62607015e-25

julia> convert_unit(1, :meV, :mK)
11604.518121550082

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

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

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

Equivalent to numpy meshgrid.