This package does quantum operator algebra (i.e., algebra with noncommuting operators) in Julia, supporting bosonic, fermionic, and twolevel system operators, with arbitrary names and indices, as well as sums over any of the indices. It defines an opinionated canonical form (normal ordering plus some additional rules) to automatically simplify expressions. It is recommended to use an interface that can display LaTeX formulas (e.g., Jupyter notebooks) for convenient output formatting.
Starting from v1.4, QuantumAlgebra also interoperates with computer algebra systems (CAS) such as Symbolics.jl or SymPy.jl / SymPyPythonCall.jl, as the "scalar" prefactors of each quantum term can be arbitrary expressions provided by these systems. While such expressions do not support symbolic indices in the same way as QuantumAlgebra, they provide much more flexibility in terms of the mathematical operations and powerful manipulation functions possible on the parameters.
Example jupyter notebooks are available in the examples
folder and can be
viewed online with
nbviewer
and tried out interactively with
Binder.
Please see the release notes for a summary of changes in each version.
The basic functions to create QuantumAlgebra expressions (which are of type
QuExpr
) are

a(inds...)
anda'(inds...)
for a and a^{†}, the annihilation and creation operators for a bosonic mode. 
f(inds...)
andf'(inds...)
for f and f^{†}, the annihilation and creation operators for a fermionic mode. 
σx(inds...)
,σy(inds...)
,σz(inds...)
for the Pauli matrices σ^{x,y,z} for a twolevel system (TLS). 
σp(inds...)
,σm(inds...)
for excitation and deexcitation operators σ^{±} for a twolevel system (TLS). 
Indices: All of these functions take an arbitrary number of indices as arguments, which can be either integers (1,2,...) or symbolic, where symbolic indices must be a single unicode character, with possibly an integer subindex:
julia> using QuantumAlgebra julia> a() a() julia> a'(:i) a†(i) julia> f'(1,2,:i_9) f†(12i₉) julia> σx(:i_1, 1, :j, :k_2, :μ_2, :◔_1, :😄_121) σˣ(i₁1jk₂μ₂◔₁😄₁₂₁)

You can define your own bosonic/fermionic/twolevel system operators with a set of macros:
@boson_ops name
defines new function$name()
(and deprecated$(name)dag()
) for bosonic speciesname
.@fermion_ops name
defines new function$name()
(and deprecated$(name)dag()
) for fermionic speciesname
.@tlsxyz_ops name
defines new functions$(name)x()
,$(name)y()
and$(name)z()
for the Pauli matrices for twolevel system speciesname
.@tlspm_ops name
defines new functions$(name)p()
and$(name)m()
for the twolevel system excitation and deexcitation operators for speciesname
.
Note that for
@boson_ops
and@fermion_ops
, the deprecated$(name)dag()
functions are defined for backward compatibility. These will be removed in a future version, as$(name)'()
is now the preferred syntax for creating an adjoint.julia> @boson_ops b (b (QuExpr constructor), b† (QuExpr constructor)) julia> b'(:k)*b(:i) b†(k) b(i)
Operators with different names are assumed to belong to different "species" and always commute. For fermions, this is not always desired, since you might want to use different named operators to refer to different kinds of states for the same species (e.g., localized and itinerant electrons). This can be achieved with the macro
@anticommuting_fermion_group
, which creates several fermionic operators that mutually anticommute:julia> @anticommuting_fermion_group c d julia> normal_form(c()*d() + d()*c()) 0

param(name::Symbol,state='n',inds...)
to create a named parameter.state
must be one of'r'
,'n'
, or'c'
for purely real, nonconjugated complex, and conjugated complex parameters. More conveniently, parameters can be entered with string macrosPr"name_inds..."
andPc"name_inds..."
for real and complex parameters:julia> Pr"g_i,j_2,k" g(ij₂k) julia> Pr"g_i,j_2,k" == param(:g,'r',:i,:j_2,:k) true julia> Pc"α_3" == param(:α,3) true

Arithmetic operations (
*
,+
,
,^
,adjoint
='
) are supported (exponents must be nonnegative integers), with anyNumber
types integrating automatically. Division by numbers is also supported.julia> 5*a'(:k)*f(3)*σx(3) 5 a†(k) f(3) σˣ(3) julia> (5//3+4im) * a'(:k)*f(3)*σx(3) + 9.4 9.4 + (5//3+4i) a†(k) f(3) σˣ(3) julia> (a(:i)*f(:k))' f†(k) a†(i)
If you explicitly need a bare number as a QuantumAlgebra expression, you can use, e.g.
QuExpr(1)
(which is equal toone(QuExpr)
). However, most functions that take aQuExpr
will also accept a bare number. 
∑(ind,A::QuExpr)
to represent an analytic sum over indexind
. Since summed indices have no semantic meaning, the index within the expression gets replaced by a special numbered sum index#ᵢ
, withi=1,2,...
.julia> ∑(:i,a(:i)) ∑₁ a(#₁)

normal_form(A::QuExpr)
converts an expression to a welldefined "canonical" order. To achieve this canonical form, relevant commutators etc are used, so an expression written as a single product can turn into a sum of expressions. The order is essentially normal ordering (creation before annihilation operators, with σˣʸᶻ in the middle), with some additional conventions to make the normal form (hopefully) unique. In some contexts (e.g., interactive work), it can be convenient to automatically transform all expressions to normal form. This can be enabled by callingQuantumAlgebra.auto_normal_form(true)
. To make the setting permanent, callQuantumAlgebra.auto_normal_form(true; set_preference=true)
or alternatively use Preferences.jl directly, i.e., callPreferences.set_preferences!(QuantumAlgebra,"auto_normal_form"=>true/false)
.julia> normal_form(a(:i)*a'(:j)) δ(ij) + a†(j) a(i)

expval(A::QuExpr)
to represent an expectation value.julia> expval(a'(:j)*a(:i)) ⟨a†(j) a(i)⟩

expval_as_corrs(A::QuExpr)
to represent an expectation value through its correlators, i.e., a cumulant expansion.julia> expval_as_corrs(a'(:j)*a(:i)) ⟨a†(j)⟩c ⟨a(i)⟩c + ⟨a†(j) a(i)⟩c

comm(A::QuExpr,B::QuExpr)
to calculate the commutator [A,B] = AB  BA.julia> comm(a(),a'()) a†() a() + a() a†() julia> normal_form(comm(a(),a'())) 1

Avac(A)
andvacA(A)
simplify operators by assuming they are applied to the vacuum from the left or right, respectively. To be precise,Avac(A)
returns A' such that A0⟩ = A'0⟩, whilevacA(A)
does the same for ⟨0A. These functions automatically applynormal_form
to assure that the operators are simplified as much as possible. Note that "vacuum" for twolevel systems is interpreted as the lower state,σᶻ0⟩ = 0⟩
.julia> Avac(a()) 0 julia> Avac(a(:i)*a'(:j)) δ(ij) julia> Avac(a()*a'()*a'()) 2 a†() julia> vacA(a()*a'()*a'()) 0 julia> Avac(σx()) σˣ() julia> Avac(σz()) 1
Both functions can also be called with an optional second argument,
Avac(A,modes_in_vacuum)
orvacA(A,modes_in_vacuum)
, which is an iterable over operators (or a single operator) that will be assumed to be in the vacuum state, while all others are not. Note that the operators inmodes_in_vacuum
do not distinguish by index, i.e., if the modes have indices, all modes with the same name are assumed to be in the vacuum state. To avoid confusion, themodes_in_vacuum
argument thus does not accept operators with indices.julia> Avac(a(),a()) 0 julia> Avac(a(),f()) a() julia> Avac(a(:i)*a'(:j),f()) δ(ij) + a†(j) a(i) julia> Avac(a'()*a()*f()*f'(),f()) a†() a() julia> @boson_ops b julia> Avac(a'()*a()*b()*b'()^2*f()*f'(),(f(),b())) 2 a†() b†() a()

vacExpVal(A,S=1)
calculates the vacuum expectation value ⟨0S^{†}AS0⟩, i.e., the expectation value ⟨ψAψ⟩ for the state defined by ψ⟩=S0⟩. The result is guaranteed to not contain any operators.julia> vacExpVal(a'()*a()) 0 julia> vacExpVal(a'()*a(), a'()^4/sqrt(factorial(4))) 4.000000000000001 julia> vacExpVal(a'()*a(), a'()^4/sqrt(factorial(big(4)))) 4 julia> vacExpVal(σx()) 0
Like
vacA
andAvac
,vacExpVal
also takes an optionalmodes_in_vacuum
argument,vacExpVal(A,S,modes_in_vacuum)
(since all arguments are positional,S
has to be given explicitly in this case even if it is just the identity operator, i.e.,vacExpVal(A,1,a())
):julia> @boson_ops b julia> vacExpVal(a'()*a()*b()^2*b'()^2*f()*f'(), 1, (f(),b())) 2 a†() a()

heisenberg_eom(A,H,Ls=())
calculates the Heisenberg equation of motion for operatorA
under the action of HamiltonianH
and potential Lindblad decay termsLs
, given by dA/dt = i[H,A] + ∑_{i} γ_{i} (L_{i}^{†} A L_{i}  1/2 {L_{i}^{†} L_{i},A}). The Lindblad decay operators are passed as a tuple (not an array) of tuples, where each inner tuple describes one decay operator. The possible forms are(L,)
for decay operatorL
,(γ,L)
for decay operatorL
with rateγ
, and(inds,γ,L)
for decay operators summed over the given indices (note that this is different from the operator itself being a sum, seen in the example below). Finally,L
can (in all three cases above) be just a single operator or a tuple of two operatorsL=(X,Y)
to represent offdiagonal Lindblad terms L_{X,Y}[ρ] = X ρ Y^{†}  1/2 {Y^{†} X,ρ}.julia> H = Pr"ω"*a'()a() julia> Ls = ((Pr"γ",a()),) julia> heisenberg_eom(a(),H,Ls) 1//2 γ a()  1i ω a() julia> H = QuExpr() julia> Ls = ((:i,a(:i)),) julia> heisenberg_eom(a(:i),H,Ls) 1//2 a(i) julia> Ls = ((∑(:i,a(:i)),),) julia> heisenberg_eom(a(:i),H,Ls) 1//2 ∑₁ a(#₁) julia> Ls = (((:i,:j),(a(:i),a(:j))),) julia> heisenberg_eom(a(:i),H,Ls) 1//2 ∑₁ a(#₁)

heisenberg_eom_system(H,rhsfilt,Ls=(),ops=nothing)
calculates the system of equations of motion for the expectation values of operators appearing inH
andLs
(same conventions as forheisenberg_eom
above). Typically, these equation systems are not closed without approximations as equations for products of n operators involve products of m>n operators, so the system has to be truncated. This is achieved with a filter function that removes higherorder terms or rewrites them (approximately) in terms of lowerorder expressions. The functionrhsfilt
is applied to the righthand side of the equations to filter them as desired. Ifrhsfilt(A::QuExpr)::QuExpr
is a function, it will be applied to the calculated righthand side of the equations.QuantumAlgebra
comes with two predefined constructors for filter functions,droplen(maxorder::Int)
, which leads to all terms of order higher thanmaxorder
being neglected, anddropcorr(maxorder::Int)
, where all terms of order higher thanmaxorder
are rewritten in terms of lowerorder expressions up to ordermaxorder
and higherorder correlators, with those correlations being neglected (i.e.,dropcorr(1)
will replace ⟨a^{†} a⟩ = ⟨a^{†} a⟩_{c} + ⟨a^{†}⟩ ⟨a⟩ ≈ ⟨a^{†}⟩ ⟨a⟩). Ifrhsfilt
is a number, it will be interpreted asdroplen(rhsfilt)
. Finally, theops
argument can be used to specify the operators that should be used to "seed" the system of equations, otherwise all operators appearing inH
are used.julia> H = Pr"ω"*a'()*a() + Pr"χ"*a'()*(a'()+a())*a(); julia> Ls = ((Pr"γ",a()),); julia> heisenberg_eom_system(H,2,Ls,a()) dₜ⟨a()⟩ = 1//2 γ ⟨a()⟩  1i ω ⟨a()⟩  2i χ ⟨a†() a()⟩  1i χ ⟨a()²⟩ dₜ⟨a†() a()⟩ = γ ⟨a†() a()⟩ dₜ⟨a()²⟩ = 2i χ ⟨a()⟩  γ ⟨a()²⟩  2i ω ⟨a()²⟩
The
heisenberg_eom_system
function can also be passed eitherExpVal
orCorr
as a first argument, which will give the equations of motion of the expectation values (the default) or correlators (corresponding to a cumulant expansion) of the operators.julia> H = Pr"ω"*a'()*a() + Pr"χ"*a'()*(a'()+a())*a(); julia> Ls = ((Pr"γ",a()),); julia> heisenberg_eom_system(Corr,H,1,Ls,a()) dₜ⟨a()⟩c = 1//2 γ ⟨a()⟩c  1i ω ⟨a()⟩c  2i χ ⟨a†()⟩c ⟨a()⟩c  1i χ ⟨a()⟩c²

julia_expression(A)
to obtain a julia expression that can be used to automatically build codes implementing equations derived with QuantumAlgebra. Every expectation value or correlator is treated as a separate array. Daggers are represented asᴴ
, which are valid identifiers that can appear in the array names. Note that expectation values and correlators are not distinguished, so it is best to have all expressions use the same kind.julia> julia_expression(expval_as_corrs(a'(:j)*a(:i))) :(aᴴ[j] * a[i] + aᴴa[j, i])
Also note that expressions are always treated as arrays, even if they have no indices (which gives zerodimensional arrays). If you are working with scalar quantities exclusively, it might be useful to clean up the resulting expression (e.g., use
MacroTools
to remove the[]
).julia> julia_expression(expval(a'()*a()*σx())) :(aᴴaσˣ[])

By default, twolevel system operators are represented by the Pauli matrices
σˣʸᶻ
, and callingσp()
andσm()
will give results expressed through them:julia> σp() 1//2 σˣ() + 1//2i σʸ() julia> σm() 1//2 σˣ()  1//2i σʸ()
This can be changed by calling
QuantumAlgebra.use_σpm(true; set_preference=true/false)
(where the value ofset_preference
determines whether this is stored permanently using Preferences.jl). In this mode,σ⁺
andσ⁻
are the "fundamental" operators, and all expressions are written in terms of them. Note that mixing conventions within the same expression is not supported, so it is suggested to set this flag once at the beginning of any calculation.julia> QuantumAlgebra.use_σpm(true) julia> σp() σ⁺() julia> σx() σ⁺() + σ⁻() julia> σz() 1 + 2 σ⁺() σ⁻()
Several preferences changing the behavior of QuantumAlgebra can be set permanently (this uses Preferences.jl):
"define_default_ops"
: if this is set tofalse
(default istrue
), the "default" operatorsa, adag, f, fdag, σx, σy, σz, σp, σm
are not defined upon import. Note that changing this value requires restarting the Julia session to take effect. The setting can be changed withQuantumAlgebra.set_define_default_ops(true/false)
(which will inform you whether a restart is required) or withPreferences.set_preferences!(QuantumAlgebra,"define_default_ops"=>true/false).
"auto_normal_form"
: Choose whether all expressions are automatically converted to normal form upon creation. The default isfalse
. It can be changed for a single session withQuantumAlgebra.auto_normal_form(true/false)
, and can be made permanent withQuantumAlgebra.auto_normal_form(true/false; set_preference=true)
or withPreferences.set_preferences!(QuantumAlgebra,"auto_normal_form"=>true/false)
. Note that this could previously be set by defining an environment variable"QUANTUMALGEBRA_AUTO_NORMAL_FORM"
, but this usage has been deprecated and will be removed in a future version."use_σpm"
: Choose whether for twolevel systems, the "basic" operators are excitation/deexcitation operatorsσ⁺
,σ⁻
or the Pauli matricesσˣ
,σʸ
,σᶻ
. This can be changed in a single session by callingQuantumAlgebra.use_σpm(true/false)
, and can be made permanent withQuantumAlgebra.use_σpm(true/false; set_preference=true)
or withPreferences.set_preferences!(QuantumAlgebra,"use_σpm"=>true/false)
.
If you use QuantumAlgebra in academic work, we would appreciate a citation. See
CITATION.bib
for the relevant references.