Operators

# Operators

Operators can be defined as linear maps from one Hilbert space to another. However, equivalently to states, operators in QuantumOptics.jl are interpreted as coefficients of an abstract operator in respect to one or more generally two, possibly distinct Bases. For a certain choice of bases $\{|u_i\rangle\}_i$ and $\{|v_j\rangle\}_j$ an abstract operator $A$ has the coefficients $A_{ij}$ which are connected by the relation

$A = \sum_{ij} A_{ij} | u_i \rangle \langle v_j |$

For this reason all operators define a left hand as well as a right hand basis:

mutable struct MyOperator{BL<:Basis,BR<:Basis} <: AbstractOperator{BL,BR}
basis_l::BL
basis_r::BR
# ...
end

For performance reasons there are several different implementations of operators in QuantumOptics.jl, all inheriting from the abstract AbstractOperator type:

They have the same interface and can in most cases be used interchangeably, e.g. they can be combined using arithmetic functions *, /, +, -:

b = SpinBasis(1//2)
sx = sigmax(b)
sy = sigmay(b)
sx + sy
sx * sy

Additionally, the following functions are implemented for all types of operators, if possible:

Conversion from one type of operator to another is also provided. I.e. to obtain a DenseOperator or SparseOperator use the DenseOperator constructor and sparse function, respectively.

### Operator data and tensor products

The data field of an operator (or a ket/bra) built by a tensor product exhibits reverse ordering to the standard Kronecker product, i.e. tensor(A, B).data = kron(B.data, A.data). This is due to the fact that this order respects the column-major order of stored data in the Julia language which is beneficial for performance. One has to keep this in mind when manipulating the data fields. If desired you can change the data output printed to the REPL with the QuantumOptics.set_printing function, i.e. by doing QuantumOptics.set_printing(standard_order=true). Note, that this will only change the displayed output while leaving the respective operator data fields the unmodified.

## Dense operators

DenseOperator is the default type used for density operators. I.e. creating an operator by using the tensor product of a ket and a bra state results in a DenseOperator. It is implemented as:

type DenseOperator{BL<:Basis,BR<:Basis,T<:Matrix{ComplexF64}} <: AbstractOperator{BL,BR}
basis_l::BL
basis_r::BR
data::T
end

where the data is stored as complex (dense) matrix in the data field.

The DenseOperator(::AbstractOperator) constructor can be used to convert other types of operators to dense operators.

## Sparse operators

SparseOperator is the default type used in QuantumOptics.jl. The reason is that in many quantum systems the Hamiltonians and jump operators in respect to the commonly used bases are sparse. They are implemented as:

using SparseArrays # hide
type SparseOperator{BL<:Basis,BR<:Basis,T<:SparseMatrixCSC{ComplexF64,Int}} <: AbstractOperator{BL,BR}
basis_l::BL
basis_r::BR
data::T
end

To convert other operators to sparse operators the sparse(::AbstractOperator) function can be used.

## Lazy operators

Lazy operators allow delayed evaluation of certain operations. This is useful when combining two operators is numerically expensive but separate multiplication with states is relatively cheap. A nice example is the transform implemented for particles. It allows using a fast fourier transformation to convert a state from real space to momentum space, applying a diagonal operator and converting it back. Doing this in operator notation is only fast if the the order of operations is $\mathrm{IFFT}*(D*(\mathrm{FFT}*\psi))$. To create a Hamiltonian that uses this calculation order, lazy evaluation is needed:

xmin = -5
xmax = 5
Npoints = 100
b_position = PositionBasis(xmin, xmax, Npoints)
b_momentum = MomentumBasis(b_position)

p = momentum(b_momentum)
x = position(b_position)

Tpx = transform(b_momentum, b_position);
Txp = dagger(Tpx)

H_kin = LazyProduct(Txp, p^2/2, Tpx)
H_pot = x^2
H = LazySum(H_kin, H_pot)

In this case the Hamiltonian $H$ is a lazy sum of the kinetic term $p^2/2$ and the potential term $x^2$ where the kinetic term is the lazy product mentioned before. In the end this results in a speed up from $O(N^2)$ to $O(N \log N)$.

Besides the above LazyProduct, there is also an implementation for lazy sums and lazy tensor products. While a LazySum works very much identical to the LazyProduct, a LazyTensor is slightly different in terms of implementation. As a brief example, consider the case of two spin-1/2 particles:

b0 = SpinBasis(1//2)
b = tensor(b0, b0)

sm0 = sigmam(b0) # Single spin operator

# Build composite space using lazy tensors
sm1 = LazyTensor(b, [1], [sm0])
sm2 = LazyTensor(b, [2], [sm0])

H = LazySum(LazyProduct(dagger(sm1), sm1), LazyProduct(dagger(sm2), sm2))

Note

A LazyTensor can only consist of SparseOperator and/or DenseOperator when it is to be used with a time evolution. Using, for example, LazyProduct to build a LazyTensor will result in an error. However, in almost all use cases, one can rewrite these constructs such that LazyTensor remains at the lowest level.