using QuantumAlgebra
# convience function to show both the original and normal_form version of an expression
dispnormal(x) = display("text/latex","\$ $(latex(x)) \\quad\\to\\quad $(latex(normal_form(x))) \$");
Create parameters with param(:name,state,indices...)
, where state = 'r','n','c'
stands for real, non-conjugated, or conjugated parameters. For more convenience, we also have the string macros Pr"name"
for real parameters and Pc"name"
for complex (non-conjugated) parameters. These macros support indices with Pr"name_i,j"
.
ωij = param(:ω,'r',:i,:j)
α = param(:α,'n')
display(ωij)
display(α)
display(ωij == Pr"ω_i,j")
display(ωij')
display(α')
true
normal_form
also reorders parameters
x = α*α*ωij*α*ωij*α'
display(x)
normal_form(x)
Parameters where state is not set are complex by default.
x = param(:g,2)*param(:g,1)*param(:g,'c',3)*param(:b)*param(:a)*param(:a)*param(:d,'c',3)*param(:f,'r',1)
dispnormal(x)
dispnormal(x')
Three types of operators are currently supported: Bosonic, fermionic, and two-level system operators. By default, QuantumAlgebra
defines a
as a bosonic operator, f
as a fermionic operator, and σx
,σy
,σz
,σp
,σm
as two-level system operators (σp
/σm
are $\sigma^\pm$). These are functions that create QuExpr
objects. Note that the adjoint '
can be applied directly to these functions, or to the resulting objects, i.e., a'() = a()'
. While the recommended interface is to use a'()
and f'()
for constructing creation operators, for backwards compatibility, the functions adag = a'
and fdag = f'
dispnormal(a()*a'())
with indices
dispnormal(a(:i)*a'(:j))
Indices can be symbolic or integer:
dispnormal(a(:i)*a'(1))
Fermions anticommute
dispnormal(f(:i)*f'(:j))
dispnormal(f(:i,:j)*f'(:i,:l))
$f^2$ and $f^{\dagger 2}$ are normalized to $0$.
dispnormal(f()^2)
dispnormal(f'()^2)
Since normal form also gives canonical ordering of operators by indices, anticommutation of fermionic operators can introduce a minus sign.
dispnormal(f(:j)*f(:i))
Indices in $\delta$s are automatically ordered and simplified (so $\delta_{ik}\delta_{kj} \to \delta_{ij}\delta_{ik}$).
dispnormal(f(:i,:k)*f'(:k,:j))
Bosonic and fermionic operators commute (independent of indices).
dispnormal(a(:i)*f(:i)*f'(:j))
Define a second bosonic species with name b
. These always commute with other species.
@boson_ops b;
dispnormal(a(:i)*b(:i)*b'(:j))
Define a second fermionic species with name g
. These always commute with any other species.
@fermion_ops g;
dispnormal(f(:i)*g(:j) + g(:j)*f(:i))
To create (mutually anticommuting) fermionic operators referring to different kinds of states of the
same species, use @anticommuting_fermion_group
.
@anticommuting_fermion_group c d e;
anticomm(A,B) = A*B + B*A;
dispnormal(anticomm(c(),d()))
dispnormal(anticomm(c(),e'()))
dispnormal(anticomm(c(:i),c'(:j)))
dispnormal(anticomm(d(:i),e'(:j)))
By default, operators for two-level systems (Pauli matrices) are written in terms of the Hermitian matrices $\sigma^x,\sigma^y,\sigma^z$ (we use superscripts so as to avoid possible confusion with indices).
display(σx(3))
display(σz(3))
display(σp(1))
normal_form
contracts products of Pauli matrices
dispnormal(σx()^2)
... but not if indices are different (since the expression would be $\delta_{ij} + \sigma^x_i \sigma^x_j (1-\delta_{ij})$)
dispnormal(σx(:i)*σx(:j))
Commutation is performed, using $[\sigma^a, \sigma^b] = 2 i \varepsilon_{a b c} \sigma^c$
dispnormal(σz(:i)*σx(:j))
Excitation/deexcitation operators $\sigma^\pm$ get rewritten in terms of $\sigma^{x,y,z}$ by default.
dispnormal(σm(:j)*σp(:i))
With QuantumAlgebra.use_σpm()
, we can ask to use the excitation/deexcitation operators $\sigma^+$ and $\sigma^-$ instead of the Pauli matrices $\sigma^{x,y,z}$ as the "natural" basis.
Note that exchanging them easily produces "complicated" expressions since they are neither bosons nor fermions, but most properly thought of as "hard-core bosons".
QuantumAlgebra.use_σpm()
dispnormal(σm(:j)*σp(:i))
Use QuantumAlgebra.use_σxyz()
to switch back to $\sigma^{x,y,z}$.
QuantumAlgebra.use_σxyz()
dispnormal(σm(:j)*σp(:i))
Scalars are incorporated "naturally", with the code trying to use integers and rationals where possible.
dispnormal(a()*a'() - 1)
dispnormal(1//2*a()*a'()*a'())
dispnormal(5//4*σx()*σy()*σz())
dispnormal(6//3*σx()*σy()*σz())
dispnormal(σp(1) * σm(1) - 1//2)
Define sums over indices with ∑(ind
, where ∑
is entered as \sum<tab>
(note that this is not Σ = \Sigma<tab>
). Sum indices are automatically renamed to the special symbol #ᵢ
(with $i=1,2,\ldots$) and reordered since they have no definite identity.
Sums over $\delta$s disappear automatically.
dispnormal(∑(:i, a(:i)*a'(:j)))
dispnormal(∑(:i, Pr"g_i"*a(:i)*a'(:j)))
H = ∑(:i,Pr"ω_i"*a'(:i,:i))
display(H)
dispnormal(a(:k,:k)*H)
H = ∑(:i,Pr"ω_i"*a'(:i)*a(:i))
display(H)
dispnormal(a(:k)*H)
x = ∑(:i,a'(:i)*a(:i))
dispnormal(x)
dispnormal(a'(:n)*x)
dispnormal(a(:n)*x)
dispnormal(σz(:n)*x)
dispnormal(Pc"g_n"*a(:n)*x)
dispnormal(x*a'(:n))
dispnormal(x*a'(:n)*a'(:n))
dispnormal(x*a'(:n)*a(:n))
dispnormal(comm(σx(5),σy(3)))
dispnormal(comm(σx(5),σx(5)))
dispnormal(comm(σx(1),σz(1)))
dispnormal(2*comm(σx(5),σy(5)))
dispnormal(comm(σy(5),comm(σx(5),σy(5))))
dispnormal(comm(2//5*param(:h)*σx(5),3*param(:g)*σy(5)))
dispnormal(comm(σp(1),σp(1)))
dispnormal(comm(σp(1),σm(1)))
dispnormal(comm(σm(1),σp(1)))
dispnormal(comm(σm(1),σm(1)))
dispnormal(comm(σm(1),σp(2)))
dispnormal(comm(σm(1),σz(1)))
dispnormal(comm(a(2),a'(1)))
dispnormal(comm(a(1),a'(1)*a(1)))
dispnormal(comm(a(1),a'(1)*a(2)*a(1)))
dispnormal(comm(a(1),a(2)*a'(1)*a(1)*a'(2)))
dispnormal(σp(1)*σz(1) + σp(1))
H = ∑(:i,Pr"ω_i"*a'(:i)*a(:i)) + ∑(:j,1//2*Pr"ωe_j"*σz(:j)) + ∑(:i,∑(:j,Pr"g_i,j"*(a'(:i)+a(:i))*σx(:j)))
tmp1 = normal_form(comm(a'(:n)*a(:m),H))
tmp2 = normal_form(comm(a'(:n)*a(:n),H))
display(tmp1)
display(tmp2)
display(tmp2 - QuantumAlgebra.replace_inds(QuantumAlgebra.QuIndex(:m)=>QuantumAlgebra.QuIndex(:n))(tmp1))
for op in [a(:n),a'(:n),σx(:k),σy(:k),σz(:k)]
display("text/latex",string("\$[",latex(op),",H] = ",latex(normal_form(comm(op,H))),"\$"))
end
for op in [a'(:n)*σz(:k),a(:n)*σz(:k)]
display("text/latex",string("\$[",latex(op),",H] = ",latex(normal_form(comm(op,H))),"\$"))
end
expval_as_corrs
expresses expectation values of operator products in terms of correlators.
expval(H)
expval_as_corrs(H)
display(expval_as_corrs(a(2)))
display(expval_as_corrs(3*a(2)))
using Latexify
display(expval_as_corrs(a(2)*a(2)))
display(expval_as_corrs(2*a(1)*a(2)*a(3)))
display(expval_as_corrs(a(1)*a(2)*a(3)*a(4)))
display(expval_as_corrs(a'(2)*a(1)*σz(1)))
display(expval_as_corrs( Pr"g"*a'(3)*a(2)*a(2)))
display(expval_as_corrs(-Pr"g_1"*σz(1)))
display(vacA(Avac(σp(1)*σm(1))))
display(vacA(Avac(σm(1)*σp(1))))
nphotstate(n,ind) = a'(ind)^n / sqrt(factorial(n))
for i2 = (:n,:m)
stateop = 1/√2*nphotstate(3,:n) + 1/√2*nphotstate(1,i2)
display("text/latex",string("\$|\\psi\\rangle = \\left(",latex(stateop),"\\right)|0\\rangle\$"))
for A in [one(QuExpr),a'(:n)*a(:n),a'(:n)*a'(:n)*a(:n)*a(:n)]
display("text/latex",string("\$\\langle\\psi|",latex(A),"|\\psi\\rangle = ",latex(vacExpVal(A,stateop)),"\$"))
end
i2 == :n && display("text/latex","")
end
For converting to code, we cannot have bare operators, but only expectation values or correlators/cumulants. An expectation value $\langle a^\dagger_i a_j\rangle$ is represented as a two-dimensional array aᴴa[i,j]
, etc. Sums are not written explicitly, but indicated by special sum index names s̄ᵢ
, which are not possible "normal" index symbols (since "s̄"
is a two-character unicode string, and we only allow single-character indices). $\delta$ are converted to I
(code using this should have using LinearAlgebra
, which defines I
as the UniformScaling
type, with I[i,j] = δᵢⱼ
).
x = ∑(:i, a(:i)*a'(:j)) + f(:k_1)*f'(:k_2)
x = expval(normal_form(x))
display(x)
julia_expression(x)
:(1.0 + I[k₁, k₂] + aᴴa[j, s̄₁] + -1.0 * fᴴf[k₂, k₁])