Hi everyone!
I’d like to announce the release of QuantumCumulants.jl, a package that automizes the symbolic derivation and numerical solution of mean-field equations.
Why it’s useful
Generalized mean-field equations are often used to overcome the size-limit of operators represented by matrices in a quantum system. For example, the Hilbert space for N two-level atoms scales in size as 2^N, so storing a state vector or operator that acts on this Hilbert space for even a few tens of atoms is already challenging.
Mean-field equations are derived by representing operators as noncommutative variables, which follow fixed fundamental commutation relations. You can use those commutation relations to derive equations of motion (cf. Heisenberg equation). When forming the average over them, you obtain a set of c-number ordinary differential equations. However, without further approximations this set is likely to be extremely large or even infinite.
To overcome this issue you can truncate the system by neglecting higher-order quantum correlations using the generalized cumulant expansion method. Using this method you can factorize average values of products of operators at an arbitrary order into lower ones, e.g. in lowest order you have \langle A B\rangle \approx \langle A \rangle\langle B \rangle for two operators A and B.
While the approach is generally quite useful, it can quickly become cumbersome and error-prone as the number of equations needed to describe a complex system at higher order orders is quite large. For example, I personally spent two weeks of my PhD studies deriving a set of 20 such equations, where the majority of the time was spent on fixing errors in the equations and the code that implemented them.
QuantumCumulants.jl does all of this for you. This makes it especially useful to treat large quantum systems with a low degree of “quantumness”. It also enables the treatment of systems at higher orders were the number of equations would already be too large to deal with analytically. Finally, it makes the numerical implementation simple by allowing you to convert the resulting set of equations to ModelingToolkit.
How it works
The framework implements its own custom noncommutative algebra. Fundamental commutation relations are hard-coded into the package and applied immediately when constructing symbolic expression. This makes it efficient and ensures correctness when later performing the cumulant expansion. The cumulant expansion is performed automatically to any order specified by the user, such that it is straightforward to obtain a closed set of equations.
Simplifications other than the application of commutation relations are done using the Symbolics.jl package. Parameters and averages are symbolic numbers and are directly implemented as Sym types from the SymbolicUtils package. That makes it very easy to convert resulting equations of motion to an ODESystem from ModelingToolkit giving you access to all its cool features and fast numerical solutions with DifferentialEquations.jl.
At this point I’d like to thank the JuliaSymbolics team, especially @shashi and @ChrisRackauckas, for helping me getting things to work with SymbolicUtils, Symbolics and ModelingToolkit!
Short example
Here’s a short example deriving and solving the first-order mean-field equations of a single-atom laser:
using QuantumCumulants
h_cav = FockSpace(:cavity)
h_atom = NLevelSpace(:atom, (:g,:e))
h = tensor(h_cav, h_atom)
@cnumbers Δ g κ γ ν
@qnumbers a::Destroy(h) σ::Transition(h)
H = Δ*a'*a + g*(a'*σ(:g,:e) + a*σ(:e,:g))
J = [a,σ(:g,:e),σ(:e,:g)]
rates = [κ,γ,ν]
eqs = meanfield([a,σ(:g,:e),σ(:e,:e)], H, J; rates=rates, order=1)
using ModelingToolkit, OrdinaryDiffEq
sys = ODESystem(eqs)
p0 = (Δ=>0, g=>1.5, κ=>1, γ=>0.25, ν=>4)
u0 = ComplexF64[1e-2, 0, 0]
prob = ODEProblem(sys,u0,(0.0,50.0),p0)
sol = solve(prob,RK4())
For more details and examples, check out the documentation.