ABSTRACT (original article): MulAtoLEG (Multi-Atom Liouville Equation Generator) is an open-source Mathematica package for generating Liouville superoperators and Liouville equations, specialized for multilevel atomic systems comprising an arbitrary number of atoms. This scheme is based on an extension to multilevel atomic systems, originally developed by Lehmberg [R. H. Lehmberg, Phys. Rev. A 2, 883 (1970)] as an adjoint master equation for ensembles of two-level emitters and later reformulated by Genes [M. Reitz, C. Sommer and C. Genes, PRX Quantum 3, 010201 (2022)] as a master equation. The package facilitates the generation of equations for complex transition configurations in alkali atoms. Although primarily designed for atomic systems, it can also generate the master and adjoint master equations for general Hamiltonians and Lindbladians. In addition, it includes functionalities to construct the differential equations in the dressed-state basis, where, in many cases, the non-unitary evolution operator can be determined explicitly. To maximize computational efficiency, the package leverages Mathematica’s vectorization and sparse linear algebra capabilities. Since MulAtoLEG produces exact equations without approximations, the feasible system size is naturally limited by the available computational resources. CITATION (original article): Pablo Yanes-Thomas, Rocío Jáuregui-Renaud Santiago F. Caballero-Benítez, Daniel Sahagún Sánchez, Alejandro Kunold (2025), MultiAtomLiouvilleEquationGenerator: A Mathematica package for Liouville superoperators and master equations of multilevel atomic systems, arXiv:2512.23591. https://doi.org/10.48550/arXiv.2512.23591​
GitHub: https://github.com/alexkunold/MultiAtomLiouvilleEquationGenerator

Content


Loading MultiAtomLiouvilleEquationGenerator

To load the MultiAtomLiouvilleEquationGenerator, write:
SetDirectory[PathToMultiAtomLiouvilleEquationGenerator];
Needs[“MultiAtomLiouvilleEquationGenerator`”]
Set the path to the MultiAtomLiouvilleEquationGenerator package and load
In[]:=
PathToMultiAtomLiouvilleEquationGenerator="/Users/akb/Dropbox/Trabajo/Rubidio/Project008";
In[]:=
SetDirectory[PathToMultiAtomLiouvilleEquationGenerator];​​Needs["MultiAtomLiouvilleEquationGenerator`"]
Content

Functions

These are the functions

MultiAtomBasis

MultiAtomBasis[na, nl] : na, nl ∈ Integers
Returns a SparseArray containing a multi atom basis for na atoms with nl quantum levels each.
Content

Example

Generate a matrix basis for two atoms with two quantum levels each
In[]:=
h=MultiAtomBasis[2,2]
Out[]=
SparseArray
Specified elements: 36
Dimensions: {16,4,4}

Show that the basis elements are orthonormal
In[]:=
Table[Tr[el1.el2],{el1,h},{el2,h}]//MatrixForm
Out[]//MatrixForm=
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
1
Generate a matrix basis for two atoms with four quantum levels
In[]:=
h=MultiAtomBasis[2,4]
Out[]=
SparseArray
Specified elements: 1296
Dimensions: {256,16,16}


Sigma

Sigma[tr,α, nl, na] : tr={
tr
1
,
tr
2
};
tr
1
,
tr
2
,α,na,nl∈
Ingeters
Returns a SparseArray with the jump operator 
tr
1
〉
〈
tr
2
 between states 
tr
1
〉
and 
tr
2
〉
for the α-th atom of a system consisting of na atoms with nl quantum levels.
Content

Example

Generate a the jump operator 
1〉
〈
2
 (t r = {1,2}) for the first atom (α=1) in a two level (nl=2) three atom (na=3) system
In[]:=
σ=Sigma[{1,2},1,2,3]
Out[]=
SparseArray
Specified elements: 4
Dimensions: {8,8}

The corresponding matrix is given by
In[]:=
σ//MatrixForm
Out[]//MatrixForm=
0
0
0
0
1
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0

SigmaT

SigmaT[tr,α, nl, na] : tr={
tr
1
,
tr
2
};
tr
1
,
tr
2
,α,na,nl∈
Ingeters
Returns a SparseArray with the adjoint jump operator 
tr
2
〉
〈
tr
1
 between states 
tr
1
〉
and 
tr
2
〉
for the α-th atom of a system consisting of na atoms with nl quantum levels.
Content

Example


Sigmaz

Sigmaz[tr,α, nl, na] : tr={
tr
1
,
tr
2
};
tr
1
,
tr
2
,α,na,nl∈
Ingeters
Returns a SparseArray with the
Pauli
σ
z
-like
operator 
tr
1
〉
〈
tr
1
 - 
tr
2
〉
〈
tr
2
 of states 
tr
1
〉
and 
tr
2
〉
for the α-th atom of a system consisting of na atoms with nl quantum levels.
Content

Example


Proj

Proj[h, Op]: h, Op ∈ SparseArray
Returns a SparseArray of dimension Length[h] = Length[Op] that contains the expansion coefficients
Op
i
of the operator Op given by
Op
i
= Tr[
h
i
Op]
where h is the matrix basis.
Content

Examples


ExpectationValue

ExpectationValue[h, q, Op]: h, q, Op ∈ SparseArrayReturns an expression for the expectation value
∑
i
Tr[
h
i
Op]
q
i
of the operator Op where
q
i
=Tr[
h
i
ρ] are the density matrix coefficients and h is the matrix basis.
Content

Example


Ham

Ham[transCE, transdd, nl, na, t] : transCE, trasdd ∈ List; na, nl ∈ Ingeters, t is a Mathematica symbol that represents the time variableReturns a SparseArray for the time-dependent Hamiltonian in the interaction picture that accounts for light-matter and dipole-dipole interactions. The variable transCE = {{{
ql
1,1
,
ql
1,2
},
α
1
,
Rb
1
,
Δ
1
,
ϕ
1
}, {{
ql
2,1
,
ql
2,2
},
α
2
,
Rb
2
,
Δ
2
,
ϕ
2
}, …} contains the parameters for the light-matter interaction of the
α
k
-th atom. These parameters are: the two quantum levels excited by the laser
ql
k,1
and
ql
k,2
, the Rabi parameter
Rb
k
, the detuning
Δ
k
and the phase
ϕ
k
. The variable transdd={{
tr
1,1
,
tr
1,2
},
α
1
,
β
1
,
Ω
1
,
δ
1
}, {
tr
2,1
,
tr
2,2
},
α
2
,
β
2
,
Ω
2
,
δ
2
}, …} contains the parameters for the dipole-dipole interaction.
tr
k,1
={
ql
k,1,1
,
ql
k,1,2
},
tr
k,1
={
ql
k,2,1
,
ql
k,2,2
} are the two interacting dipole transitions corresponding to the two distinct
α
k
-th and
β
k
-th atoms (
α
k
≠
β
k
). The strength of the interaction is given by
Ω
k
.
δ
k
is the detuning between these two transitions. For
δ
k
to be significant, it must be much smaller than both energy transitions themselves. Usually
δ
k
=0. transCE and transdd might be empty lists in case there are no light-matter interactions or/and there are no dipole-dipole interactions.
Content

Example


HamDressed

HamDressed[transCE, transdd, nl, na] : transCE, trasdd ∈ List; na, nl ∈ IngetersReturns a SparseArray for the time-independent Hamiltonian in the dressed state base that accounts for light-matter and dipole-dipole interactions. The variable transCE = {{{
ql
1,1
,
ql
1,2
},
α
1
,
Rb
1
,
Δ
1
,
ϕ
1
}, {{
ql
2,1
,
ql
2,2
},
α
2
,
Rb
2
,
Δ
2
,
ϕ
2
}, …} contains the parameters for the light-matter interaction of the
α
k
-th atom. These parameters are: the two quantum levels excited by the laser
ql
k,1
and
ql
k,2
, the Rabi parameter
Rb
k
, the detuning
Δ
k
and the phase
ϕ
k
. The variable transdd={{
tr
1,1
,
tr
1,2
},
α
1
,
β
1
,
Ω
1
,
δ
1
}, {
tr
2,1
,
tr
2,2
},
α
2
,
β
2
,
Ω
2
,
δ
2
}, …} contains the parameters for the dipole-dipole interaction.
tr
k,1
={
ql
k,1,1
,
ql
k,1,2
},
tr
k,1
={
ql
k,2,1
,
ql
k,2,2
} are the two interacting dipole transitions corresponding to the two distinct
α
k
-th and
β
k
-th atoms (
α
k
≠
β
k
). The strength of the interaction is given by
Ω
k
.
δ
k
is the detuning between these two transitions. For HamDressed these parameters must be
δ
k
=0, otherwise the Hamiltonian becomes time-dependent. transCE and transdd might be empty lists in case there are no light-matter interactions or/and there are no dipole-dipole interactions.

RotatingFrameParms

RotatingFrameTrans

LiouvilleLindbladian

LiouvilleLindbladianDagger

LiouvilleCommutator

LiouvilleCommutator[h, H] : h, H ∈ SparseArray
Returns a SparseArray of dimension Length[h] that contains the Liouville superoperator corresponding to the commutator i[H, ρ] of the master equation or to the commutator i[H, Q] of the adjoint master equation. It is important to notice that in the case of the master equation it is important to set the sign of the Liouville superoperator to negative. In contrast, in the adjoint master equation, this adjustment is not necessary.

Example

Assuming that the transitions 1〉 → 2〉 and 2〉 → 3〉 have different energy differences and, therefore, do not interact among themselves, the transdd parameter is given by
The Hamiltonian is generated with
From this Hamiltonian ham we generate the Liouville superoperator of the commutator
Note the negative sign at the beginning of the instruction.

LiouvilleMasterEquation

LiouvilleMasterEquation[q, Li, t] : q, Li ∈ SparseArray, t is a Mathematica symbol that represents the time variable
Returns a List of dimension Length[q] = Length[Li] that contains the Liouville differential equations that correspond to the Liouville superoperator Li which includes both the Lindbladian and Hamiltonian parts generated using LiouvilleCommutator and LiouvilleLindbladian. The variable q should be a list of the density matrix coefficients, i.e., Tr[h ρ(t)] .

LiouvilleAdjointMasterEquation

MasterEquation

AdjointMasterEquation

CleanUp

CleanUp[Op] : Op ∈ SparseArray
Returns a SparseArray of dimension Length[Op] with a cleaner set of ArrayRules. For example, in some cases it might eliminate vanishing imaginary elements.

FilteredClebschGordan

DipoleMatrixElement

DipoleMatrixElement[i, jp, j, fp, f, mfp, mf ] i, jp, j, fp, f, mfp, mf ∈ Integers
Returns the angular part of dipole matrix element 〈i, jp, fp, mfp p  i, j, f, mf 〉 leaving out the reduced matrix element. This separation follows the Wigner–Eckart theorem, in which a matrix element is expressed as the product of a reduced matrix element (independent of magnetic quantum numbers) and a Clebsch–Gordan coefficient or 3j-symbol that carries the angular dependence. i and j are the nuclear and electronic angular momentum quantum numbers, respectively, while f and mf denote the hyperfine total angular momentum and its magnetic quantum number, respectively.

InitialConditions

LiouvilleLindbladianFromLindbladOperators

LiouvilleDaggerLindbladianFromLindbladOperators

TransitionLists

FarField

WDipoleDipoleElement

MDipoleDipoleElement

Examples

These are the examples

A single two-level atom

We set up the master equation for a single two-level atom. Using its solution, we plot the time dependence of the quantum-level populations. We also compute the emission as a function of the laser detuning.

Two two-level atoms

We set up the master equation for two two-level atoms. Using its solution, we plot the time dependence of the quantum-level populations. We also perform a spectroscopic analysis of the emission by sweeping the detuning and evaluating the correlation between the emissions of the two atoms.

Five two-level atoms

We set up the master equation for a chain of five two-level atoms coupled through nearest-neighbor interactions. Using its solution, we plot the time dependence of the correlations between the emissions of neighboring atoms and of the atoms at the edges of the chain. We also perform a spectroscopic analysis of the emission by sweeping the detuning and evaluating these correlations.

Optical pumping

Example

To distinguish the dynamics of all the populations, we repeat the previous example by slightly shifting the polarization vectors
Load MultiAtomLiouvilleEquationGenerator
We set the parameters and build the master equation (see previous example)

Superradiance

We illustrate the phenomenon of superradiance by comparing the spontaneous emission of a single atom with that of five coupled atoms. Superradiance is a collective quantum effect where a group of emitters radiate light much more intensely and faster than they would if each one emitted independently.

The dressed state evolution operator

We illustrate how to transform the solution of the master equation from the dressed-state basis to the standard-state basis for a system of two three-level atoms. This method provides a way to obtain the explicit form of the evolution operator. The resulting solution is shown to be consistent with that obtained directly from the master equation in the standard basis.

Example

Load MultiAtomLiouvilleEquationGenerator
This is an important step because, although the evolution operator can be obtained symbolically for very sparse Liouville operators, in most cases it must be calculated numerically. The transCE configuration variable, which describes this light–matter interaction, is
Assuming that the transitions 1〉 → 2〉 and 2〉 → 3〉 have distinct energy gaps and therefore do not interact with each other, the transdd parameter is given by
Calculate the dressed state Hamiltonian and the corresponding Liouville operator
Calculate the Lindbladian superoperator
It is important to remember that the Lindbladian superoperator remains invariant under rotating-frame transformations and, therefore, is calculated just as it would be in the standard basis. The full Liouville operator is the sum of LiouvilleHR and LiouvilleL
Obtain the rotating-frame transformations
and use them to convert the results from the rotating frame to the standard basis
To illustrate the procedure for calculating expectation values, we construct the operator op = 1〉〈2⊗1+ 2〉〈1⊗1, and then evaluate its expectation value
To compare, we also find the system’s evolution through the master equation
and numerically solve the differential equations
and find the expectation value of op using this solution
Superimpose the solutions obtained from the dressed state and the standard base

Time evolution of transmons

To illustrate the use of LiouvilleLindbladianFromLindbladOperators and LiouvilleDaggerLindbladianFromLindbladOperators in obtaining the dynamics of an arbitrary system, we analyze the time evolution of a transmon subject to both a coherent drive and dissipative mechanisms.

Example

Load MultiAtomLiouvilleEquationGenerator
where ns is the number of states. We analyze the behaviour of the first 8 states of a transmon. We thus start by generating a matrix basis with 8 states with
We also define the ascent and descent operators for such a number of states and a SparseArray called one for the unitary matrix
The Hamiltonian is thus given by
We set the constants of the system
and calculate the eigenvalues of H in order to know the excitation energy Δ between the ground and the first excited states
We calculate the Liouville operator for the Hamiltonian
The Lindbladian Liouville operator is calculated from just one Lindblad operator, namely a, and just one decay coefficient γ
Compute the overall Liouville operator and build the differential equations
Build the system of differential equations for the density matrix coefficients
We reset the consts variable to include the decay coefficient $\gamma$ and a drive with frequency Δ that excites the transition from the ground state to the first excited state. We then evaluate the equations with these updated constants and set the initial conditions such that the transmon starts in the ground state. Finally, the initial conditions are combined with the equations.
Solve the differential equations
Calculate and plot the populations of the 8 transmon quantum levels
We observe that as the ground state (blue) is gradually depopulated, the first excited state (orange) becomes populated due to the driving field. In addition, a small fraction of the population leaks into the higher excited states (green and subsequent levels), reflecting the multi-level structure of the transmon and the fact that the drive is not perfectly restricted to a two-level subspace.

CITE THIS NOTEBOOK

Multi atom Liouville equation generator​
by Pablo Yanes-Thomas, Rocío Jáuregui-Renaud Santiago F. Caballero-Benítez, Daniel Sahagún Sánchez & Alejandro Kunold​
Wolfram Community, STAFF PICKS, January 27, 2026
​https://community.wolfram.com/groups/-/m/t/3629205