PRX QUANTUM 6, 010201 (2025)

Tutorial

Solving Quantum Dynamics with a Lie-Algebra Decoupling Method

1,2,* 2,3,†
Sofia Qvarfort and Igor Pikovski
1
Nordita, KTH Royal Institute of Technology and Stockholm University, Hannes Alfvéns väg 12,
SE-106 91 Stockholm, Sweden
2
Department of Physics, Stockholm University, AlbaNova University Center, SE-106 91 Stockholm, Sweden
3
Department of Physics, Stevens Institute of Technology, Castle Point on the Hudson, Hoboken,
New Jersey 07030, USA

(Received 12 September 2023; revised 27 September 2024; published 27 January 2025)

Quantum technologies rely on the control of quantum systems at the level of individual quanta. Math-
ematically, this control is described by Hamiltonian or Liouvillian evolution, requiring the application of
various techniques to solve the resulting dynamic equations. Here, we present a tutorial for how the quan-
tum dynamics of systems can be solved using a Lie-algebra decoupling method. The approach involves
identifying a Lie algebra that governs the dynamics of the system, enabling the derivation of differential
equations to solve the Schrödinger equation. As background, we include an overview of Lie groups and
Lie algebras aimed at a general-physicist audience. We then prove the Lie-algebra decoupling theorem
and apply it to both closed and open dynamics. The results represent a broad methodology to find the
dynamics of quantum systems with applications across many fields of modern quantum research.

DOI: 10.1103/PRXQuantum.6.010201

CONTENTS V. NONLINEAR QUANTUM DYNAMICS:

CAVITY OPTOMECHANICS 16
I. INTRODUCTION 1 VI. DECOUPLING THEOREM FOR
II. INTRODUCTION TO LIE ALGEBRAS OPEN-SYSTEM DYNAMICS 17
AND LIE GROUPS 3 A. Introduction to vectorization 18
A. Lie groups 3 B. Vectorizing the Lindblad equation 18
B. Lie algebras 3 VII. OPEN DYNAMICS: A THERMALIZING
C. Link between Lie groups and Lie algebras 4 HARMONIC OSCILLATOR 19
III. THE LIE-ALGEBRA DECOUPLING VIII. SUMMARY 20
THEOREM 4 ACKNOWLEDGMENTS 20
A. Proof of the decoupling theorem 6 APPENDIX: PARTITIONING QUANTUM
B. Lie-algebra decoupling theorem in phase DYNAMICS 21
space 8 1. Unitary dynamics 21
C. A concise recipe for decoupling 9 2. Nonunitary dynamics 22
IV. EXAMPLES OF LINEAR (GAUSSIAN) REFERENCES 22
DYNAMICS 10
A. Hamiltonian with linear terms 10 I. INTRODUCTION
B. Hamiltonian with quadratic terms 12 Quantum physics results in unique phenomena such as
C. Most general single-mode Gaussian quantum superpositions and entanglement, which can be
Hamiltonian 15 harnessed for applications ranging from quantum tech-
nologies to tests of fundamental physics. However, devel-
oping such applications requires exact knowledge and
* control on the level of individual quanta and their time
Contact author: sofi[Link]@[Link]
†
Contact author: [Link]@[Link] evolution. As new experimental platforms and hybrid sys-
Published by the American Physical Society under the terms of
tems emerge, novel dynamics and interaction regimes
the Creative Commons Attribution 4.0 International license. Fur- become increasingly accessible. In order to fully unlock
ther distribution of this work must maintain attribution to the the potential of these systems, it is crucial to be able to
author(s) and the published article’s title, journal citation, and accurately model their quantum dynamics. To mention just
DOI. Funded by Bibsam. a few examples, quantum dynamics plays a crucial role for

2691-3399/25/6(1)/010201(24) 010201-1 Published by the American Physical Society

SOFIA QVARFORT and IGOR PIKOVSKI PRX QUANTUM 6, 010201 (2025)

quantum control [1], quantum information processing [2], solutions, the system dynamics can be solved exactly. The
and quantum sensing [3]. Beyond quantum technologies, Lie-algebra method has recently been employed to study
searches for new effects in fundamental physics frequently a number of different problems, including the nonlinear
lead to predicted changes in the quantum dynamics of a dynamics of optomechanical systems [19–25], the dynam-
system. To detect these often extremely weak effects, it is ics of coupled oscillators [26,27], cooling protocols [28],
crucial to be able to model the system dynamics exactly. time crystals in nonlinear Kerr cavities [29], Floquet engi-
As an example, the search for a quantum theory of gravity neering [30], and state preparation [31], as well as the study
has resulted in the study of modifications to the usual alge- of open-systems dynamics [24,32–34].
bras used in quantum theory [4,5]. Such a modified algebra The goal of this tutorial is to provide an introduction to
necessarily leads to changes in the dynamics of quantum the Wei-Norman Lie-algebra decoupling method and how
systems. In order to predict observable effects, we typically it can be used to solve quantum dynamics. To build intu-
require methods to construct the resulting decoupled uni- ition for the method, we provide an introduction to Lie
tary operations, at least perturbatively [6]. These methods groups and Lie algebras aimed at a general-physicist audi-
are also relevant for the study of nonclassicality [7] and ence. We then state and prove the Lie-algebra decoupling
nonlinearities [8] in quantum systems. theorem, as well as an analogue version in phase space.
It is, however, generally challenging to treat quantum To demonstrate its applicability, we apply the theorem to a
dynamics analytically. The core of the difficulty lies in the number of examples, including time-dependent linear and
noncommutativity of operators that enter into the Hamil- quadratic Hamiltonians, as well as a nonlinearly coupled
tonian. While it is always possible to express the formal bipartite system. In addition, we detail how the Lie-algebra
solution to the dynamics in terms of an exponential opera- decoupling theorem can be applied toward solving quan-
tor, this operator cannot always be tractably applied to the tum master equations, such as the Gorini-Kossakowski-
initial quantum state and, therefore, the dynamics cannot Sudarshan-Lindblad equation [35,36] (often just referred
easily be studied. The notion of solving quantum dynam- to as the Lindblad equation). We do so by vectorizing the
ics can therefore be understood as finding a closed-form state and writing the Lindblad equation as a linear matrix
expression of the time-evolution operator that facilitates equation [37]. In this paper, we focus on applications in
the computation of any quantity of the system, such as quantum optics and, in particular, on continuous-variable
expectation vales. Ultimately, this is equivalent to solv- quantum optical systems. But we highlight that the math-
ing Schrödinger’s differential equation, which can be a ematical techniques are widely relevant for other fields as
challenging task. well. To give a few brief examples in other fields, in many-
In many fields of quantum science, it is often sufficient body physics, such techniques are useful to solve classes
to consider only small quantum perturbations around clas- of bosonic models, such as particles interacting with inde-
sical solutions or to average over many systems. In such pendent bosons [38]. These models constitute some of the
cases, perturbation theory can be used to derive approxi- few exactly solvable many-body problems and thus have
mate solutions to the dynamics. As long as the effects of many applications, such as for the description of relaxation
interest are weak, these perturbative solutions often result effects [16], or electron gases and quantum liquids [39], to
in reasonably accurately models of experiments in the lab- name a few. In nuclear physics, some problems on vibra-
oratory. Indeed, many mathematical methods have been tional motion can be mapped onto bosonic Lie algebras
developed to express, manipulate, and truncate exponen- [40], such that the decoupling methods discussed here are
tial operators, including the Zassenhaus formula [9], the of relevance. Lie-algebra methods are also applied in NMR
Magnus expansion [10], exponential-expansion methods physics [41]. In addition, related methods known as uni-
[11], with a variety of results presented in Ref. [12], and tary integration [42–44] have been successfully applied to
the Suzuki-Trotter decomposition [13–16] (for a review, treat finite-dimensional systems with both closed [45,46]
see also Ref. [17]). In order to control the full quantum and open dynamics [47]. More broadly, the Lie-algebra
behavior of individual quantum systems, however, such decoupling method can become useful whenever a phys-
as in quantum optics, quantum information, and for quan- ical problem (potentially with many interacting systems)
tum technologies, it can be necessary to go beyond these can be mapped onto a relatively simple closed algebra.
approximations. The tutorial is structured as follows. Section II pro-
One mathematical technique that can be used to analyt- vides a mathematical introduction to Lie groups and Lie
ically solve quantum dynamics has been put forward by algebras. Section III contains a proof of the decoupling
Wei and Norman [18]. At its core is the observation that a theorem, as well as a straightforward recipe for how it
Lie algebra generated by a Hamiltonian can be used as a can be applied. Then, in Sec. IV, we apply the theorem to
basis for studying the dynamics. Expanding the time evo- solve the dynamics of a Hamiltonian with time-dependent
lution of the system in terms of this basis then allows us linear and quadratic interaction terms, respectively. We
to derive a set of scalar differential equations. In the cases also consider an example of nonlinear dynamics in Sec. V
in which the resulting differential equations have analytic (where by nonlinear, we refer to Hamiltonian terms of

010201-2
SOLVING QUANTUM DYNAMICS. . . PRX QUANTUM 6, 010201 (2025)

creation and annihilation operators with powers larger Let us now focus specifically on Lie groups. They
than 2), namely, the dynamics that arise from the non- are continuous groups and they play a ubiquitous role
linear optomechanical Hamiltonian. Next, we consider the in physics and mathematics. For example, in quantum
application of the Lie-algebra decoupling theorem to open- mechanics, the set of unitary time-evolution operators
system dynamics in Sec. VI, as well as an example of a form a Lie group, as we will see below. We proceed with
thermalizing harmonic oscillator in Sec. VII. the formal definition of a Lie group.

Definition (Lie group). A Lie group is a set G with two

II. INTRODUCTION TO LIE ALGEBRAS AND LIE
structures:
GROUPS
The Lie-algebra decoupling theorem can be understood (1) G is a group with the structure discussed in the
and applied without deep knowledge of Lie algebras but, definition of a group shown above.
for the interested reader, here we provide a basic mathe- (2) G is a smooth and real manifold. Smoothness means
matical introduction to Lie groups and Lie algebras, as well that the group operation and inverse map are dif-
as the connection between them. In each case, we start the ferentiable. The group is described by a set of real
discussion with a formal definition and then provide a few parameters that describe the group elements.
common examples.
Examples of Lie groups include, e.g., the translation group,
the special unitary group SU(n), the group of all invert-
A. Lie groups ible linear maps, and the special orthogonal group in three
Groups arise in many separate contexts in physics. To dimensions, SO(3). For all of these groups, the inverse
understand Lie algebras, it helps to first introduce the con- and the zero elements can be constructed. The real line
cept of Lie groups. We start by recounting the general has continuous elements, where the inverse can be con-
formal definition of a group. structed by subtraction, and the identity element is zero.
The special unitary group SU(n) describes n × n unitary
matrices with determinant 1. Here, the inverse element can
Definition (Group). A group (G, ∗) is a set G with group be obtained through complex conjugation and the iden-
elements g and a binary operation ∗ such that G × G → G, tity is the n × n identity matrix. Finally, the SO(3) group
which satisfies three conditions: describes the rotation of vectors in three dimensions. The
inverse element can be constructed through orthogonal-
(1) Associativity. For any three elements x, y, z ∈ G, we
ity and matrix multiplication automatically satisfies the
have (x ∗ y) ∗ z = x ∗ (y ∗ z).
requirement of associativity. In fact, most continuous rota-
(2) Identity. The group must include an identity element
tion groups are Lie groups. A list of Lie groups can be
 ∈ G. It must hold that ∀g ∈ G, multiplying g by
found in Ref. [48].
the identity element leaves g invariant, such that  ∗
g = g ∗  = g.
(3) Inverse. For each element in the group, there must B. Lie algebras
be an inverse element. That is, for each element g ∈ We are now finally ready to properly define a Lie alge-
G, there is some inverse element ḡ ∈ G such that bra. Lie algebras are the key object of interest in this
g ∗ ḡ = ḡ ∗ g = . work, and the link between Lie groups and Lie algebras
underpins the Lie-algebra decoupling theorem. Given a Lie
Perhaps one of the simplest examples of a group is the set group, it is always possible to construct a Lie algebra from
of discrete rotations that leave a square invariant. That is, the group and sometimes there are advantages to studying
the group elements correspond to rotation actions (acting the algebra rather than the group itself. We begin with the
on a square) and the group operation corresponds to the basic definition of a Lie algebra [49].
composition of actions. This group is known as Z4 and
contains rotations of zero, 90◦ , 180◦ , and 270◦ . This set Definition (Lie algebra). A Lie algebra is a vector space
of rotations satisfies the group axioms, since any combi- g over some field [50] F, together with a binary operation
nation of the group elements satisfies the same symmetry. [·, ·] : g × g → g (the Lie bracket), which must satisfy the
That is, we can rotate the square by first 90◦ , then 180◦ , and following axioms:
it remains invariant. The associativity criterion is satisfied
by a proper representation of the group elements, such as a (1) bilinearity, such that [ax + by, z] = a[x, z] + b[y, z]
matrix. Lastly, the identity is equivalent to not rotating the and [z, ax + by] = a[z, x] + b[z, y] for all scalars
square and the inverse can easily be constructed by rotat- a, b in the field F and all elements x, y, z in g
ing the square in the reverse direction by the same number (2) alternativity, which means that the Lie bracket is
of degrees. zero for the same element: [x, x] = 0

010201-3
SOFIA QVARFORT and IGOR PIKOVSKI PRX QUANTUM 6, 010201 (2025)

(3) the Jacobi identity, which states that III. THE LIE-ALGEBRA DECOUPLING
THEOREM
[x, [y, z]] + [z, [x, y]] + [y, [z, x]] = 0 (1)
Equipped with some knowledge of Lie groups and Lie
In fact, the commutator bracket [A, B] = AB − BA, which algebras, we are now ready to study the Lie-algebra decou-
is commonly used in quantum physics, satisfies these pling theorem, originally outlined by Wei and Norman
criteria. This will be important to us later. [18]. This section closely follows the proof developed in
Ref. [18] but with slightly different notation in order to
C. Link between Lie groups and Lie algebras be consistent with modern conventions in quantum infor-
The next question is how Lie groups connect with Lie mation and quantum optics—for a presentation of these
algebras. While this can be discussed in great mathemat- methods in the context of optomechanical systems, see also
ical detail, here we provide an example that is hopefully Ref. [51]. For convenience, we set  = 1 in this section.
intuitive to the quantum physicist. In short, the Lie alge- Intuitively, the Lie-algebra decoupling method can be
bra generates the group. To begin with, let us consider a thought of as separating a dynamical problem into the
Lie group L with elements G(α) ∈ L, where α is some real notion of directions of evolution and the strength of the
parameter. To determine the action of the element near the evolution, where the directions are defined by the alge-
identity, we can slightly perturb G(α) for a small αj , to find bra elements and the strength corresponds to how strongly
each algebra element acts on the quantum state. For exam-
G(α) ≈ 1 + iδαj Xj , (2) ple, in continuous-variable quantum information, rotation,
displacement, and squeezing operators are often used to
where we have defined the generator Xj . Then, perform- describe the trajectory of a quantum state in phase space.
ing this small perturbation many times in addition to the A state can, e.g., be strongly squeezed but only slightly
identity operation, we find displaced away from the vacuum state. We differentiate
  between the squeezing and displacement operators and the
iαj Xj k extent to which they are applied.
lim 1 + ≡ eiαj Xj = D(α), (3)
k→∞ k The Lie-algebra decoupling method effectively trans-
lates the problem of solving an operator-valued linear
where D(α) is defined as the displacement. We can now differential equation into that of solving a coupled sys-
also define the generator Xj as the rate of change with tem of differential equations. One advantage of using this
respect to the parameter αj : method is that problems that would have required the use
 of large numerical Hilbert spaces can instead be treated
∂ 
Xj ≡ −i D(α) . (4) by solving a set of (potentially coupled) scalar differential
∂αj j equations. Crucially, it also allows us to solve the dynam-
ics of time-dependent Hamiltonians, which can otherwise
We note that Eq. (3) is, in fact, the definition of the expo-
be quite tricky. While the resulting scalar differential equa-
nential map. The Xj are generators, which form a Lie
tions do not always have analytic solutions and might sim-
algebra. The Lie algebra then generates the group together
ilarly have to be solved using numerical methods, errors
with the real parameters αj . That is, given a Lie alge-
due to the limited size of numerical Hilbert spaces can be
bra with a set of n elements, one can always use the
avoided.
exponential map to generate a Lie group.
In short, the Lie-algebra decoupling method is con-
Here, we also introduce the notion of a Lie-algebra
cerned with solving the Schrödinger equation. Consider
basis. The basis is made up of the Lie-algebra generators.
the first-order differential equation
For example, the Pauli matrices σ̂x , σ̂y , and σ̂y form a basis
together with the identity element 1 that can be used to
represent qubits and their evolution. dÛ(t)
= −iĤ (t)Û(t), (5)
Finally, before moving on, let us also touch on the notion dt
of a dimension of a Lie algebra. The Lie algebra forms
a linear vector space and the dimension of a Lie algebra where the Hermitian operator Ĥ (t) is the Hamiltonian
refers to the dimension of this vector space. However, a Lie (which we have assumed is time dependent, to make it as
algebra with finite dimension can have an infinite number general as possible) and Û(t) is a time-evolution operator.
of elements, since any arbitrary linear combination of the The formal solution to Û(t) is given by
elements also belongs to the Lie algebra. This should be   
contrasted with the dimension of the Hilbert space, which i t
Û(t) = T exp − dt Ĥ (t) , (6)
is the representation space. To summarize, the evolution  0
of a quantum system that lives in an infinite-dimensional
Hilbert space could be generated by a finite Lie algebra. where T indicates time ordering.

010201-4
SOLVING QUANTUM DYNAMICS. . . PRX QUANTUM 6, 010201 (2025)

We then assume that the Hamiltonian Ĥ (t) can be writ- where â and ↠are annihilation and creation operators that
ten as a finite sum with m terms. A Hamiltonian with satisfy the canonical commutation relation [â, ↠] = 1. The
an infinite number of unique terms would by extension term (↠â)2 is often referred to as a Kerr nonlinearity. The
also generate an infinite-dimensional Lie algebra. The Lie- dynamics generated by Ĥ can be trivially solved since all
algebra decoupling theorem holds for finite-dimensional terms commute. However, if we now add a linear drive
Lie algebras. We therefore assume that the initial Hamil- term, such that Ĥ = ω↠â + g(↠â)2 + ξ ↠+ ξ ∗ â, where
tonian can be written as a sum over m constant operators ξ is the drive-strength amplitude, we generate infinitely
Ĥ j multiplied by time-dependent coefficients Gj (t): many unique terms when we commute the Kerr nonlinear-
ity with the drive term. As a result, the system dynamics

m cannot be solved exactly and we must resort to pertur-
Ĥ (t) = Gj (t)Ĥ j . (7) bation theory. A similar situation arises in the context of
j =1 modified quantum theories, such as in the study of the gen-
eralized uncertainty principle [4,5], and the dynamics in
The set {Ĥ j } with j = 1, 2, . . . , m reproduces the Hamil- such models are often treated perturbatively [6]. Neverthe-
tonian Ĥ (t). It can be extended to a larger set with n ≥ m less, the methods presented here still remain useful as long
elements by taking the commutator of the elements in {Ĥ j } as truncations can be justified.
and adding the result to the set of Hamiltonian terms. We We now show that the existence of a finite Lie algebra L
can then write the original Hamiltonian as enables the decoupling of the time-evolution operator Û(t)
into a product of n operators, namely,

n
Ĥ (t) = Gj (t)Ĥ j , (8) Û(t) = Û1 (t)Û2 (t) . . . Ûn (t), (10)
j =1
where each component operator Ûj (t) is an operator satis-
where the coefficients with j > m are set to zero so that fying
we still retain the original Hamiltonian. Extending the sum
in this way makes it easier to write down some relations d
further on. Note that this decomposition of the Hamiltonian Ûj (t) = −iḞj Ĥ j Ûj (t) (11)
dt
and identification of the terms Ĥ j is not necessarily unique.
We may, e.g., choose either Hermitian or non-Hermitian and where the functions Fj are functions that we wish to
Ĥ j terms. determine.
The full set of Hamiltonian terms {Ĥ 1 , Ĥ 2 , Ĥ 3 , . . . , Ĥ n } In quantum theory, the advantage of writing Û(t) in the
form a finite n-dimensional Lie algebra L under commu- form given in Eq. (10) is that when the action of each Ûj (t)
tation. In other words, it is possible to find the full Lie is known, it becomes straightforward to apply them to a
algebra generated by Ĥ (t) by commuting all the elements quantum state in the Schrödinger picture, or an operator
Ĥ i of Ĥ (t). The Lie bracket in this case is the commuta- in the Heisenberg picture (by, e.g., inserting completeness
tor relation [Ĥ i , Ĥ j ] ≡ Ĥ i Ĥ j − Ĥ j Ĥ i . The Lie algebra L relations in between the unitaries). The advantage of such a
is constructed from all operators in Ĥ (t), plus all the Lie method over numerical solvers that use finite-dimensional
products matrices is significant, as the key task shifts from evolving
the operator-valued Û(t) to obtaining analytic expressions
for the scalar Fj functions.
Ĥ α1 , Ĥ α2 , Ĥ α3 , . . . Ĥ αr−1 , Ĥ αr · · · , (9) We are now at a point at which we can concisely state
the Lie-algebra decoupling theorem.
where αi = 1 to m, plus all linear combinations of such
products. Theorem 1 (Lie-algebra decoupling theorem). Suppose
If, through consecutive commutation, we find a finite that the linear operator Ĥ (t) can be expressed in the form
number of elements, the Lie algebra has a finite number
of elements. Such a finite algebra can always be found if 
m
the Ĥ j are finite-dimensional matrices. However, there are Ĥ (t) = Gj (t)Ĥ j , (12)
cases in which the Lie algebra is infinite, for which com- j =1
mutation of two or more operators continuously produces
new elements that are not already part of the algebra. In where m is a finite integer and where the functions Gj (t)
those cases, the dynamics can rarely be solved exactly and are scalar functions of time t and the Ĥ j are time-
perturbation theory must instead be used. Consider, e.g., a independent and unique operators defined in a Hilbert
bosonic system with a Hamiltonian Ĥ = ω↠â + g(↠â)2 , space H. Note that the operators Ĥ j do not have to be
where ω is a frequency and g is a coupling constant, and themselves Hermitian, even if we require that the full

010201-5
SOFIA QVARFORT and IGOR PIKOVSKI PRX QUANTUM 6, 010201 (2025)

Hamiltonian Ĥ (t) is Hermitian. Furthermore, the dimen- where Ŷ ∈ L. Then, we define powers of this equation as
sion of H can be either finite or infinite. Let the Lie algebra the nested operators
L generated by Ĥ (t) be of finite dimension n. Then, there
exists a neighborhood of t = 0 in which the solution of the (adX̂ )2 Ŷ = [X̂ , [X̂ , Ŷ]], (17)
equation
and so on. Thus the BCH formula can be stated as
dÛ(t)
= −iĤ (t)Û(t), (13)
dt eX̂ Ŷe−X̂ = (eadX̂ )Ŷ. (18)

with the initial condition Û(0) = 1, may be expressed in Proof of Lemma 1. We begin by defining a function
the form
∞
1
Û(t) = exp −iF1 (t)Ĥ 1 exp −iF2 (t)Ĥ 2 . . . F̂(a) = eaX̂ Ŷe−aX̂ = Ĉn an , (19)
n=0
n!
× exp −iFn (t)Ĥ n (t) , (14)
where the Ĉn are operator coefficients, which are indepen-
where Ĥ 1 , Ĥ 2 , . . . , Ĥ n is a basis for L and the set {Fj (t)} dent of a. When a = 1, the coefficients correspond to the
are scalar functions of time t. The functions Fj (t) depend case we are considering. Our goal is to derive expressions
only on the Lie algebra L and the initial functions Gj (t). for these coefficients in the form of a recursion relation.
The coefficients Fj are complex in general. However, if We first note that
the basis for the operators Ĥ j is Hermitian, then the coef- d
ficients Fj (t) must be real to ensure unitarity. There is, F̂(a) = X̂ , F̂(a) . (20)
however, no need a priori to choose a Hermitian basis da
for Ĥ j . Certain problems could, e.g., benefit from choos- Inserting Eq. (19) into Eq. (20), we find
ing the annihilation and creation operators â and ↠as the
basis. These facts have been previously highlighted in the ∞
  1 ∞
1
literature [44]. Ĉn an−1 = X , Ĉn an . (21)
n=1
(n − 1)! n=0
n!
The same decoupling of an evolution operator can also
be performed for the symplectic matrices when consid- The sum on the left-hand side can be rewritten by setting
ering the evolution under any quadratic Hamiltonian. We n → n + 1, such that we find the formula
demonstrate this fact in Sec. III B.
Ĉn+1 = X̂ , Ĉn an . (22)
A. Proof of the decoupling theorem
Our goal is to prove the decoupling theorem and show In this way, all coefficients can be generated through
how the functions Fj (t) are calculated. The proof is repeated commutation with X . We also have that C0 = Y,
based on two lemmas: the first is the well-known Baker- which follows from simply Taylor expanding the exponen-
Campbell-Hausdorff (BCH) lemma and the second one tials in Eq. (19). The coefficients in Eq. (22) can then be
concerns the closure of the Lie algebra. We begin with used to generate all the coefficients in Eq. (15). 
Lemma 1, which states the following.
The last lines in the lemma follow from the definition of
Lemma 1 (Baker-Campbell-Hausdorff). If two opera-
(adX ), as can be seen by writing
tors X̂ , Ŷ ∈ L, then eX̂ Ŷe−X̂ ∈ L and
 1
X̂
e Ŷe −X̂ 1
= Ŷ + [X̂ , Ŷ] + [X̂ , [X̂ , Ŷ]] eX̂ Ŷe−X̂ = eadX̂ Ŷ = (adX̂ )n Ŷ
2! n
n!
1 1
+ [X̂ , [X̂ , [X̂ , Ŷ]]] + · · · . (15) = Ŷ + [X̂ , Ŷ] + X̂ , [X̂ , Ŷ]
3! 2!
1
We define the new operator adX̂ , where adX̂ , X̂ ∈ L by the + X̂ , X̂ , [X̂ , Ŷ] + .... (23)
3!
equation
This concludes the proof of Lemma 1. We proceed with the
(adX̂ )Ŷ = [X̂ , Ŷ], (16) second lemma.

010201-6
SOLVING QUANTUM DYNAMICS. . . PRX QUANTUM 6, 010201 (2025)

Lemma 2 (Lie-algebra basis). Let Ĥ 1 , Ĥ 2 , . . . , Ĥ n be a Û(t) with respect to time t, we find the expression
basis for the Lie algebra L. Then, it follows that j −1 
dÛ(t)  n 
⎛ ⎞ ⎛ ⎞ = −i Ḟj (t) exp −iFk Ĥ k
r 
1 dt j =1 k=1
⎝ exp −iFj Ĥ j ⎠ Ĥ k ⎝ exp iFj Ĥ j ⎠ ⎛ ⎞
j =1 j =r n
× Ĥ j ⎝ exp −iFk Ĥ k ⎠ . (28)

n
= −i ξjk Ĥ j , (24) k=j

j =1
We then use the fact that dÛ(t)/dt = −iĤ (t)Û(t) (which
holds even when Û(t) requires time ordering) and multiply
where r = 1, . . . , n and where each ξjk ≡ ξjk (F1 , . . . , Fr ) is Eq. (28) by the inverse operator Û−1 (t) on the right-hand
a function of all its arguments. side and set the expression equal to Ĥ (t) in Eq. (12) to find
 j −1 
Proof of Lemma 2. To prove Eq. (24), we repeatedly 
n 
n 
apply Lemma 1 to the left-hand side. We demonstrate the Gj (t)Ĥ j = −i Ḟj (t) exp −iFk Ĥ k
first few lines of this proof. Consider Eq. (24) for r = 1, j =1 j =1 k=1
⎛ ⎞
which is the simplest case. We find, by using Eq. (15), 
1
× Ĥ j ⎝ exp iFk Ĥ k ⎠
k=j −1
exp −iF1 Ĥ 1 Ĥ k exp iF1 Ĥ 1
 j −1 

n 
(−iF1 ) 2
= −i Ḟj (t) exp −iFk adĤ k Ĥ j ,
= Ĥ k − iF1 [Ĥ 1 , Ĥ k ] + [Ĥ 1 , [Ĥ 1 , Ĥ k ]]
2! j =1 k=1
(−iF1 )3 (29)
+ [Ĥ 1 , [Ĥ 1 , [Ĥ 1 , Ĥ k ]]] + . . . . (25)
3!
where we have used the Baker-Campbell-Hausdorff lemma
(Lemma 1) in the second line.
Now, since the Lie algebra is closed under commutation,
By then applying Lemma 2 to the last line of Eq. (29),
this means that one of the terms eventually reads
we find

[Ĥ 1 , (adĤ 1 )n Ĥ k ] = [Ĥ 1 , Ĥ 1 ] = 0, (26) 

n 
n 
n
Gk (t)Ĥ k = −i Ḟj (t)ξkj (F1 , . . . , Fn )Ĥ k , (30)
k=1 j =1 k=1
for some integer n. As a result, Eq. (25) contains a finite
number of terms with different powers of F1 . Since all F1 where we recall that the ξkj are analytic functions of one or
are analytic, the resulting functions ξjk in the right-hand more of the functions Fj (t).
side of Eq. (24) are also analytic. The same argument can We now note that the operators Ĥ k are linearly indepen-
be made for multiplication of additional terms when r = 1. dent and effectively form a basis as part of the Lie algebra.
This concludes the proof of Lemma 2.  We use the linear independence of Ĥ k to find linear rela-
tions between the Gk (t) functions and the Ḟj (t). They are
related by the elements ξkj , which we can collect into a
We are now ready to prove the Lie-algebra decoupling
transformation matrix ξ . We define the vector of Hamilto-
theorem. The proof makes use of both Lemmas 1 and 2.
nian coefficients G = (G1 , G2 , . . . , Gn )T and the vector of
differentiated F functions Ḟ = (Ḟ1 , Ḟ2 , . . . , Ḟn )T , such that
Proof. (Lie-algebra decoupling theorem). We work
with the extended Hamiltonian as in Eq. (8): G = −iξ (F1 , . . . , Fn )Ḟ. (31)

Eq. (31) encodes n differential equations and explicitly

 n
Ĥ˜ (t) = Gj (t)Ĥ i , (27) reads
j =1 ⎛ ⎞ ⎛ ⎞⎛ ⎞
G1 ξ11 ξ12 · · · ξ1n Ḟ1
⎜G2 ⎟ ⎜ξ21 ξ22 · · · ξ2n ⎟ ⎜Ḟ2 ⎟
⎜
⎟
where n is the dimension of the Lie algebra, and consider ⎜ . ⎟ = −i ⎜ . .. .. ⎟ ⎜ .. ⎟ . (32)
⎝ .. ⎠ ⎝ .. .
..
. ⎠
. ⎝.⎠
the ansatz in Eq. (14), which states that Û(t) can be written
as a product of operators Ûj (t). When we then differentiate Gn ξn1 ξn2 ··· ξnn Ḟn

010201-7
SOFIA QVARFORT and IGOR PIKOVSKI PRX QUANTUM 6, 010201 (2025)

The differential equations can be decoupled if the matrix number of N modes. Our first step is to define a Hamilto-
ξ is invertible. The determinant det{ξ } is nonzero at t = 0, nian matrix H that acts as a mapping between these vectors
because ξ (0) = 1. Thus, ξ is invertible in some neighbor- of operators to the Hamiltonian operator. Any quadratic
hood of t = 0. Beyond t = 0, we know that ξ is generated Hamiltonian Ĥ can therefore be written as
by the analytic functions Fj (t) and must therefore also be
an analytic function of t. Generally, the analyticity and 1 T
Ĥ (t) = X̂ H(t)X̂. (34)
uniqueness of Fj (t) as solutions to a set of differential 2
equations are determined by the Picard-Lindelöf theorem
[52]. Recall that the Fj coefficients and, by extension, the ξ The time evolution of the operators X̂ is given by
matrix are determined by the initial choices of Hamiltonian
coefficients Ĥ j . If a different composition of the original X̂(t) = S(t)X̂, (35)
Hamiltonian Ĥ is used, this also changes the form of the
Fj coefficients and ξ . where S(t) is a symplectic matrix given by
Since ξ is invertible, we can write   t 
S(t) = T exp  dt H(t ) . (36)
Ḟ = f (G, F) = iξ −1 (F1 , . . . , Fn )G. (33) 0

The result is a set of n coupled differential equations with Here, T again indicates time ordering of the exponential,
boundary conditions F(0) = 0, which come from the fact  is the symplectic form that encodes the commutator rela-
that Û(t = 0) = 1. This ensures that the system has a tions, and H(t) is the Hamiltonian matrix in Eq. (34). In the
unique solution. However, we note that the differential basis of the annihilation and creation operators (â, ↠), the
equations themselves might not always have analytic solu- single-mode symplectic form  is given by
tions and may need to be solved numerically. While this
can also be challenging, it can be advantageous to solve 
N  
i 0
these scalar equations of motion compared with modeling = 1 , where 1 = , (37)
0 −i
an infinite-dimensional quantum systems in a numerically n=1
truncated Hilbert space.
This concludes the proof of the decoupling theorem.  where N is the number of modes under consideration. In
other bases, such as the position√ and momentum basis √
(X̂ , P̂), where X̂ = (↠+ â)/ 2 and P̂ = i(↠− â)/ 2,
B. Lie-algebra decoupling theorem in phase space the symplectic form is instead given by
In the study of quantum continuous variables, Gaussian
 
states, such as coherent states, squeezed states, and ther- 0 1
mal states, play key roles in many contexts [53]. Crucially, 1 = . (38)
1 0
Gaussian states are completely characterized by their first
and second moments, which means that we can model their The reason that we in this work choose to adopt the {â, ↠}
dynamics by focusing exclusively on their representation basis is because it is generally easier to commute the
in terms of first and second moments. For an introduction free-evolution term ↠â with the interaction terms in the
to quantum continuous variables, see Ref. [54]. Hamiltonian and thereby predict the effects of time evolu-
The Lie-algebra decoupling theorem can be applied in tion. However, both choices lead to equivalent results. The
phase space specifically to model the evolution of Gaus- link between the dynamics in the Hilbert space and phase
sian states. It is sometimes easier to solve the dynamics in space is
this way, since the problem of computing nontrivial com-
mutators and multiplications by congruence is reduced to
Û† (t)X̂Û(t) = S(t)X̂. (39)
that of matrix multiplication. Examples of such cases have,
e.g., been explored in a relativistic setting [55] and used
Just as with the Hilbert-space method, we can make an
to study squeezing and cooling in optomechanical systems
ansatz for the solution of S(t):
[21,28].
Here, we prove the Lie-algebra decoupling theorem
S(t) = S1 (t)S2 (t)S3 (t) . . . Sn (t)
for Gaussian dynamics and Gaussian states, mean-
ing that we focus exclusively on Hamiltonians with = eF1 H1 eF2 H2 eF3 H3 . . . eFn Hn , (40)
terms containing at most a quadratic number of oper-
ators. We start by defining a multimode vector of where n is the number of elements in the Lie algebra.
† † † †
first moments X̂ = (â1 , â2 , â3 , . . . âN , â1 , â2 , â3 , . . . , âN )T , Note here that our Lie algebra does not just consist of the
where â1 , â2 , â3 . . . âN are the annihilation operators for a Hamiltonian matrices Hj but, rather, the product of these

010201-8
SOLVING QUANTUM DYNAMICS. . . PRX QUANTUM 6, 010201 (2025)

matrices with the symplectic form: Hj . If we then con- Then, using the linear independence of the matrices, we
sider the individual decoupled contributions to Û(t), the find
relationship reads

n
Gk (t) = Ḟj (t)ξkj . (47)
†
Ûj (t)X̂Ûj (t) = Sj (t)X̂. (41) j =1

This means that the coefficients Fj that we have defined For systems with linear or quadratic Hamiltonian opera-
in Eq. (14) have a one-to-one relationship with those in tors, this phase-space treatment is equivalent to the gen-
Eq. (40). eral method presented in Sec. III. The main advantage is
Let us now prove the decoupling theorem in phase that it reduces a complicated operator-based multiplica-
space. We start by extending the Hamiltonian matrix to tion by congruence to that of simple matrix multiplication.
include all n elements in the algebra, such that Whether this phase-space treatment of the Hilbert-space
treatment is preferable depends on the problem at hand.

n
H(t) = Gj (t)Hj . (42) C. A concise recipe for decoupling
j =1 Here, we provide a summary of the decoupling meth-
ods in the form of a simple recipe that can be applied to
We then differentiate the ansatz in Eq. (40) with respect to any Hamiltonian that generates a finite-dimensional Lie
time t. We find algebra:
j −1  (1) Write the Hamiltonian Ĥ (t) as
dS(t)  
n
= Ḟj (t) exp[Fk Hk ]
dt j =1 k=1

m

⎛ ⎞ Ĥ (t) = Gj (t)Ĥ j (48)

n j
× Hj ⎝ exp[Fk Hk ]⎠ . (43)
k=j and identify the functions G(t) and the operators Ĥ j .
(2) Identify the Lie algebra L that is generated by
We then multiply Eq. (43) by the inverse operator S−1 (t) the Hamiltonian by commuting the Hamiltonian
on the right-hand side and use the fact that dS(t)/dt = terms Ĥ j until they produce a set of n operators
H(t)S(t) to find that is closed under commutation. For example,
[Ĥ 1 , Ĥ 2 ] ∝ Ĥ 3 . Subsequently, [Ĥ 1 , Ĥ 3 ] ∝ Ĥ 4 , and
 j −1  so on, until all algebra elements are known.

n 
n  (3) State the ansatz for the time-evolution operator,
Gj (t)Hj = Ḟj (t) exp[Fk Hk ]
j =1 j =1 k=1
⎛ ⎞ Û(t) = exp −iF1 (t)Ĥ 1 exp −iF2 (t)Ĥ 2 · · ·

1
× Hj ⎝ exp[−Fk Hk ]⎠ . (44) × exp −iFn (t)Ĥ n , (49)
k=j −1

where the Fj coefficients are dimensionless time-

By again defining dependent functions and the generators Ĥ j are part
of the Lie algebra L.
(adX )X = [X , Y], (45) (4) Differentiate the ansatz in Eq. (49) with respect to
time t to find
we can use Lemmas 1 and 2 to arrive at an expression very
d  n
similar to the one in Eq. (30): Û(t) = −iḞ1 Ĥ 1 Û(t) − iḞ2 Û1 Ĥ 2 Ûj
dt j =2

n 
n 
n
Gk (t)Hk = Ḟj (t)ξkj (F1 , . . . , Fn )Hk . (46) 
n

k=1 j =1 k=1
− iF3 Û1 Û2 Ĥ 3 Ûj
j =3

This again shows us that the algebra elements are not the 
n

bare Hamiltonian matrices Hj but, rather, the products of + · · · − iḞn Ûj Ĥ n . (50)
the Hamiltonian matrices with the symplectic form Hj . j =1

010201-9
SOFIA QVARFORT and IGOR PIKOVSKI PRX QUANTUM 6, 010201 (2025)

(5) Multiply Eq. (50) by Û−1 (t) on the right and set the Hamiltonian-interaction term. Such terms can correspond
ansatz equal to the original Hamiltonian Ĥ (t): to a number of effects. Most commonly in optical systems,
they represent continuous pumping, which arises, e.g., by
Ĥ (t) = Ḟ1 Ĥ 1 + Ḟ2 Û1 Ĥ 2 Û−1
1
injecting laser light into a cavity [56]. If the laser light
enters at a frequency different from the free frequency of
+ Ḟ3 Û2 Û1 Ĥ 3 Û−1 −1
1 Û2 + . . . . (51) the system, then the pump term changes as a function of
time. More generally, such terms are also referred to as
Evaluate all multiplications of congruence to find a drive terms.
closed-form expression. In this section, we use the Lie-algebra decoupling
(6) Use the linear independence of {Ĥ j } to construct a method to solve the dynamics of a single quantum har-
set of differential equations, where the solutions for monic oscillator with a linear driving term. The Hamilto-
Fj depend on the original Hamiltonian coefficients nian for such a system is given by
Gj (t). The equations are given by
 
Ĥ (t) = ω↠â + ω g(t)↠+ g ∗ (t)â , (53)

n
Gk (t) = −i Ḟj (t)ξkj , (52)
j =1 where ω is the free-oscillation frequency and where g(t) is
a time-dependent complex driving coefficient.
where the functions ξkj are obtained through the As a first step, we rescale time t by the frequency ω,
multiplications by congruence shown in Eq. (51). such that ωt → t, where t is now dimensionless. Then, we
(7) Solve the equations in Eq. (52) analytically or start commuting the operators in the terms of Eq. (53),
numerically for the Fj coefficients and use the result using the commutator relations [â, ↠] = 1, [↠â, â] = −â,
to determine the time-evolution operator Û(t). [↠â, ↠] = ↠. As a result, we see that the algebra that
generates the evolution of this Hamiltonian is given by
The same procedure can be carried out in phase space
(see Sec. III B). There are many Hamiltonians that can ↠â, â, ↠, 1, (54)
be treated by using this recipe. In the following sections,
we demonstrate how the dynamics of a harmonic oscilla-
where 1 is the identity operator.
tor Hamiltonian with linear and quadratic interaction terms
We then state the ansatz for the evolution operator Û(t).
can be solved.
It reads
IV. EXAMPLES OF LINEAR (GAUSSIAN) † †
DYNAMICS Û(t) = e−iFI 1 e−iF0 â â e−iF+ â e−iF− â , (55)

In this section, we study three examples of Hamiltoni-

where we have introduced the coefficients FI , F0 , F± ∈ C
ans that can be solved using the Lie-algebra decoupling
as dimensionless functions of time. Here, we note that
method. These examples are well known in quantum
since the basis of the algebra is non-Hermitian, we must
optics and have been previously considered and solved
also include the identity operator 1 in the ansatz. Given
throughout the literature. Here, we apply the Lie-algebra
a Hermitian basis, FI would be a real coefficient, which
decoupling theorem and derive the differential equations.
results in a global phase when included in the ansatz in
For full generality, we include time dependence of the
Eq. (55). Aspects of such a global phase have previously
Hamiltonians.
been discussed in Refs. [45,57].
We call the dynamics Gaussian because the Hamilto-
Before we proceed to use the Lie-algebra decoupling
nian contains at most quadratic products of annihilation
method to derive the coefficients F0 and F± , it is helpful
and creation operators. These types of dynamics are often
to note the following useful identities. First, we note that
also referred to as linear, because the resulting equations
for arbitrary complex numbers θ+ and θ− , the following
of motion for the quantum operators are linear in terms of
relations hold:
operators. Indeed, the algebra of quadratic Hamiltonians of
bosonic systems is guaranteed to have a finite Lie algebra † −θ ∗ â † ∗ 1 ∗ ∗ † 1 ∗
since the commutator of two terms with quadratic bosonic eθ+ â − = eθ+ â e−θ− â e− 2 θ+ θ− = e−θ− â eθ+ â e 2 θ+ θ− .
operators always removes one operator. (56)

A. Hamiltonian with linear terms For θ+ = θ− , this expression reduces to the familiar
† ∗
One of the simplest additions to a freely evolving relation for displacement operators D(z) = ezâ −z â =
† ∗ 2
quantum harmonic oscillator is a linear single-mode ezâ e−z â e−|z| /2 . It also directly follows from Eq. (56)

010201-10
SOLVING QUANTUM DYNAMICS. . . PRX QUANTUM 6, 010201 (2025)

that find

d
Û(t)Û−1 (t) = −iḞI − iḞ0 ↠â
† †
e−θ â âeθ â = â + θ ,
(57) dt
 
eθ â ↠e−θ â = ↠+ θ. − iḞ+ e−iF0 ↠− iḞ− eiF0 â + iF+ .
(62)
The relations in Eq. (57) are often referred to as dis-
placements of the annihilation and displacement operators, Then, we set this expression equal to −iĤ , where Ĥ is
which do not alter the commutator relation [â, ↠] = 1. the Hamiltonian with linear terms in Eq. (53). Through
In contrast, the number operator ↠â induces rotations of the linear independence of the operators, we identify the
the creation and annihilation operators. We also note that following differential equations:
any function f (â, ↠) that can be Taylor expanded changes
according to 0 = ḞI + iḞ− F+ , 1 = Ḟ0 ,
(63)
g(t) = Ḟ+ e−iF0 , g ∗ (t) = Ḟ− eiF0 .
−θ ↠â θ ↠â
e f (â, ↠)e = f (âeθ , ↠e−θ ), (58)
The second equation can be straightforwardly solved to
which holds for arbitrary θ . These quantities are useful to find F0 = t with boundary condition F0 (t = 0) = 0. This
us going forward. term corresponds to the free evolution of the harmonic
Our goal now is to determine F0 and F± in Eq. (55). We oscillator and is not changed by the linear-interaction
begin by differentiating Eq. (55) with respect to time t to terms. We then rearrange the last two equations in Eq. (63)
find to find
 t  t
d † † F+ = dt g(t )eit , F− = dt g ∗ (t )e−it . (64)
Û(t) = −iḞI e−iFI 1 e−iF0 â â e−iF+ â e−iF− â 0 0
dt
† †
− iḞ0 ↠âe−iF0 â â e−iF+ â e−iF− â From this, we note that indeed F+ = F−∗ , which is to be
† † expected given our choice of Hamiltonian. These solutions
− iḞ+ e−iF0 â â ↠e−iF+ â e−iF− â
mean that the last coefficient becomes
† †
− iḞ− e−iF0 â â e−iF+ â âe−iF− â . (59)  
t t
∗ −it
FI = −i dt g (t )e dt g(t )eit . (65)
−1 0 0
We then multiply Eq. (59) by Û (t) on the right to find
To evaluate the integrals in Eq. (64), we must first choose
d † † a specific form of the driving function g(t). Once we have
Û(t)Û−1 (t) = −iḞI − iḞ0 ↠â − iḞ+ e−iF0 â â ↠eiF0 â â done so, we can fully characterize the system dynamics.
dt
† † † † We now note that the evolution operator in Eq. (55)
− iḞ− e−iF0 â â e−iF+ â âeiF+ â eiF0 â â . (60) can be further simplified by combining the last two expo-
nentials into a displacement operator. By using the BCH
We then compute the expressions that arise from the con- formula (see Sec. III A) and the fact that F− = F+∗ , we find
gruence multiplications in the second and third terms in
† ∗ 2 /2
Eq. (60). We find, using the relations in Eqs. (57) and (58), e−iF+ â e−iF+ â = D̂(−iF+ )e|F+ | , (66)

† † â
e−iF0 â â ↠eiF0 â = e−iF0 ↠, (61a) where we have defined the displacement operator D̂(ξ ) =
† ∗
† † â eξ â −ξ â . The full evolution operator can therefore be
e−iF0 â â âeiF0 â = eiF0 â, (61b) written as
† †
e−iF+ â âeiF+ â = â + iF+ , (61c) †
−iF− â † iF− â
Û(t) = e−iϕ1 e−iF0 â â D̂(−iF+ ), (67)
e â e = â − iF− .
†
(61d)
where we have defined ϕ = FI + i|F+ |2 /2. That is, the full
We do not immediately need Eq. (61d) but it is useful to us dynamics are captured by a phase term [58], a rotation, and
later. By using Eqs. (61a), (61b), and (61c) in Eq. (60), we a displacement of the quantum state. Examining the phase

010201-11
SOFIA QVARFORT and IGOR PIKOVSKI PRX QUANTUM 6, 010201 (2025)

term, we note that it is given by For such a choice, the integrals in Eq. (64) evaluate to
 t  t
F+ =g0 [i − i cos(t) + sin(t)] ,
ϕ = −i dt dt g ∗ (t )g(t )e−i(t −t )
(75)
0 0 F− =g0 [−i + i cos(t) + sin(t)] .
 t  t
i
+ dt dt g(t )g ∗ (t )e−i(t −t )
. (68) To visualize the quadratures, we plot X̂ (t) and P̂(t)
2 0 0
parametrically as a function of rescaled time t in Fig. 1(a)
We can write the second integral as two integrals over as a function of time t for α = 1. We note that as the value
triangular regions, namely, of g0 increases, the system performs larger and larger tra-
 t  t jectories in phase space; however, it always returns to its
dt dt g(t )g ∗ (t )e−i(t −t ) initial state whenever t is a multiple of 2π .
0 0 If, instead, the function g(t) changes in time, the system
 t  t behavior becomes much more involved. Here, we find that
= dt dt g(t )g ∗ (t )e−i(t −t )
interesting effects such as resonances markedly affect the
0 0
  dynamics of the system. By resonance, we refer to time-
t t
dependent effects that occur at a frequency equal to the free
+ dt dt g(t )g ∗ (t )e−i(t −t ) . (69) frequency ω.
0 0
To explore the resonant case, we let g(t) = g0 cos(t + φ),
Then, putting the two expressions together and swapping where g0 is again the amplitude and φ is a phase off-
the labels of the second integral in Eq. (69), we find set. By solving the integrals in Eq. (64), we find that the
 t  t coefficients F± become (for φ = 0)
ϕ= dt dt Im[g ∗ (t )g(t )e−i(t −t ) ], (70)
1  1 
0 0
F+ = 1 + 2it − e2it , F− = 1 − 2it − e−2it .
4 4
which is always real and can therefore be considered a (76)
global phase.
To gain some intuition for the system dynamics, we con- We note that, compared with the coefficients in Eq. (75),
sider the evolution of the position X̂ and momentum P̂ which arise for a constant coupling, the F+ now increase
phase-space quadratures. They are given, in terms of the linearly in time.
annihilation and creation operators, by We again plot X̂ (t) and P̂(t) for the resonant driv-
ing as a function of time t. The result can be found in
↠+ â ↠− â
X̂ = √ and P̂ = i √ . (71) Figure 1(b), for different driving strengths g0 and the phase
2 2 choice φ = 0. We note that as t increases, the system is
In the Heisenberg picture, these quadratures evolve as exploring larger and larger trajectories in phase space and
X̂ (t) = Û† (t)X̂ Û(t) and P̂(t) = Û† (t)P̂ Û(t). Given an ini- does not return to its initial state.
tially coherent state |α and using the expressions in
Eqs. (61a)–(61d), we find that B. Hamiltonian with quadratic terms
We now consider harmonically trapped systems with
1     
X̂ (t) = √ eiF0 ↠+ iF− + e−iF0 â − iF+ . (72) additional quadratic Hamiltonian terms. Such terms can be
2 engineered by, e.g., changing the trapping frequency of the
Similarly, system [59]. In cases in which the quadratic terms are mod-
ulated at twice the free frequency, the term is known as a
i      parametric drive [60]. In fact, modulating the potential at
P̂(t) = √ eiF0 ↠+ iF− − e−iF0 â − iF+ . (73) parametric resonance for a specific phase offset causes a
2
reduction in the number of quanta in a harmonic oscillator
When the system starts in, e.g., a coherent state |α, the [28,60] and can in certain cases enhance the sensitivity of
expectation values are given by a quantum force sensor [23].
1     The Hamiltonian for a quantum harmonic oscillator with
X̂ (t) = √ eiF0 α ∗ + iF− + e−iF0 (α − iF+ ) , quadratic single-mode interaction terms reads
2
(74)
i     Ĥ (t) = ω↠â + ωλ+ (t)â†2 + ωλ− (t)â2 , (77)
P̂(t) = √ eiF0 α ∗ + iF− − e−iF0 (α − iF+ ) .
2
where ω is the free-oscillation frequency of the mode
Let us compute the coefficient for two specific cases. We and where the λ± (t) are complex, time-dependent, and
start by considering a constant and real coupling g(t) ≡ g0 . dimensionless coefficients.

010201-12
SOLVING QUANTUM DYNAMICS. . . PRX QUANTUM 6, 010201 (2025)

(a) (b)
6 g0 = 0.0 6 g0 = 0.0
g0 = 0.5 g0 = 0.5
g0 = 1.0 g0 = 1.0
4 4

2 2

0 0
P

P
–2 –2

–4 –4

–6 –6

–6 –4 –2 0 2 4 6 –6 –4 –2 0 2 4 6
X X

FIG. 1. Phase-space trajectories for a Hamiltonian with a linear-interaction term. Both plots show the quadratures X̂ (t) and P̂(t) as
a function of time for a constant linear term and a time-dependent linear term. When there is no driving (g± = 0), the state explores a
circle in phase space (red curve). (a) The quadrature trajectories for a constant coupling g± ≡ g0 . The state explores a limited trajectory
in phase space. (b) A time-coupling term at mechanical resonance with g± ≡ g0 cos(t + φ), which causes the state to explore wider
and wider spirals in phase space, for φ = 0 in this case. Both plots use the coherent state parameter α = 1.

To solve the dynamics induced by the Hamiltonian in a rotation with ↠â. A single-mode squeezing operator is
Eq. (77), we start by defining the following quadratic †2 ∗ 2
defined as Ŝ(ζ ) = e(ζ â −ζ â )/2 , with the complex param-
operators: eter ζ = reiϕ that includes the strength and phase-space
1 1 †  1 direction of squeezing.
K̂+ = â†2 , K̂0 = 2â â + 1 , K̂− = â2 . (78) We start by differentiating Eq. (80) with respect to time
2 4 2 t, to find
These operators form an SU(1,1) algebra and obey the
˙ = −iξ̇ K̂ e−iξ+ K̂+ e−iξ0 K̂0 e−iξ− K̂−
Û(t)
following commutation relations: + +

[K̂0 , K̂± ] = ±K̂± , [K̂+ , K̂− ] = −2K̂0 . (79) − iξ̇0 e−iξ+ K̂+ K̂0 e−iξ0 K̂0 e−iξ− K̂−

This SU(1,1) algebra shares many properties with the − iξ̇− e−iξ+ K̂+ e−iξ0 K̂0 K̂− e−iξ− K̂− . (81)
commonly used SU(2) algebra, which induces, e.g., the
dynamics of two-level systems. If SU(2) can be thought of We then multiply Eq. (81) by Û−1 (t) on the right, to find
as a sphere, the SU(1,1) algebra instead represents the two
˙ Û−1 (t) = −iξ̇ K̂ −iξ̇ e−iφ+ K̂+ K̂ eiξ+ K̂+
Û(t)
semispheres that make up the full sphere. For more back- + + 0 0
ground on the applications of SU(1,1) in quantum physics,
see Ref. [61]. See also Ref. [62] for a tutorial on the − iξ̇− e−iξ+ K̂+ e−iξ0 K̂0 K̂− eiξ0 K̂0 eiξ+ K̂+ . (82)
application of these methods to quantum optics problems.
From the last two terms in Eq. (82), we need to compute
We proceed to state the following ansatz for the time
the following multiplications by congruence:
evolution generated by the Hamiltonian in Eq. (77):
e−iξ+ K̂+ K̂0 eiξ+ K̂+ = K̂0 + iK̂+ ξ+ ,
Û(t) = e−iξ+ K̂+ e−iξ0 K̂0 e−iξ− K̂− , (80)
e−iξ0 K̂0 K̂− eiξ0 K̂0 = eiξ0 K̂− , (83)
where ξ± and ξ0 are complex time-dependent coefficients
that we wish to solve for. Note that this is one of many e−iξ+ K̂+ K̂− eiξ+ K̂+ = K̂− +2iK̂0 ξ+ −K̂+ ξ+2 .
orderings. The different orderings have been studied in
Ref. [63], where it has also been found that some of them After dividing Eq. (82) by −i and rearranging, we find
produce differential equations that correspond to the clas-
sical equations of motions of the system. We also note that ˙ Û−1 (t) = ξ̇ +iξ̇ ξ −ξ̇ ξ 2 eiξ0  K̂
Û(t) + 0 + − + +
we can reorder the exponentials in Eq. (80) into an expres-  
sion that includes a single-mode squeezing operator and + ξ̇0 + 2iξ̇− ξ+ e K̂0 + ξ̇− eiξ0 K̂− . (84)
iξ0

010201-13
SOFIA QVARFORT and IGOR PIKOVSKI PRX QUANTUM 6, 010201 (2025)

We then set this expression equal to −iĤ , where Ĥ is To simplify Eq. (89), we consider the following addition
the quadratic Hamiltonian in Eq. (77). By using the lin- formula:
ear independence of the operators, we are able to derive tan(A) + tan(B)
the following three differential equations: tan(A + B) = . (91)
1 − tan(A) tan(B)

λ+ = ξ̇+ +iξ̇0 ξ+ −ξ̇− ξ+2 eiξ0 , This, and noting that tan(iA) = i tanh(A), allows us to
write Eq. (89) as
1 = ξ̇0 + 2iξ̇− ξ+ eiξ0 , (85)  
λ− = ξ̇− eiξ0 . i 1 − 2i tanh(t)
ξ+ = −1 (92)
2λ− 1 + 21
i tanh(t)
where we recall that λ0 and λ± are dimensionless functions
We then multiply out the denominator in Eq. (92) and
of time that appear in the Hamiltonian in Eq. (77).
rearrange the expression to find
We can then rearrange the equations in Eq. (85) to
 
isolate the derivatives. We start by noting that the third 1 2 tanh(t) + 2 1
tanh(t)
equation implies that ξ̇− = λ− e−iξ0 . This allows us to ξ+ = . (93)
rewrite the first and second equations in Eq. (85) as 2λ− 1 + 2
1
i tanh(t)
 
Finally, we note that 12 + 2 2 = 12 + 2 λ+ λ− − 14 =
λ+ = ξ̇+ +iξ̇0 ξ+ −λ− ξ+2 ,
(86) 2λ+ λ− , which means that we can write Eq. (93) as
1 = ξ̇0 + 2iλ− ξ+ .  
λ+ sinh(t)
ξ+ = . (94)
We then rearrange the second equation in Eq. (86) to find  cosh(t) + 2i sinh(t)
ξ̇0 = 1 − 2iλ− ξ+ . Inserting this into the first equation in
Eq. (86) and rearranging again, we find which is the final result. The solution for ξ0 can be tested
by inserting the result in Eq. (95) into the second equation
in Eq. (85). The result satisfies the differential equation.
ξ̇+ =λ+ −iξ+ −λ− ξ+2 . (87)
The same can be done for ξ− and ξ0 , to find
 
In summary, the differential equations for ξ0 , ξ+ , and ξ− λ± sinh(t)
are given by ξ± = ,
 cosh(t) + 2i sinh(t)
  (95)
ξ̇+ = λ+ (t) − iξ+ −λ− (t)ξ+2 , i
ξ0 = −2i ln cosh() + sinh() ,
2
ξ̇0 = 2iλ− (t)ξ+ −1, (88)
ξ̇− = λ− (t)e−iξ0 , The resulting decoupling is very useful and well known in
quantum optics [65].
Our method presented here, however, directly gener-
where we have restored the potentially explicit time depen-
alizes to time-dependent coefficients. In this case, and
dence of λ± (t). Here, the equation for ξ+ is a Riccati
depending on the form of the time dependence, the differ-
equation of quadratic nonlinearity. Its nonlinearity in ξ+
ential equations in Eq. (88) must usually be solved numer-
(not to be confused with nonlinearities in terms of the
ically. There are, however, cases in which they reduce
quantum operators) is a consequence of exponentiation
to well-known differential equations, such as the Mathieu
being a nonlinear operation [64].
equation [21,28].
When the coefficients λ± (t) are constant in time with
Before moving on, we note that the quadratic Hamil-
λ± (t) = λ± , we can find an exact solution to the differen-
tonian in Eq. (77) can be cast as a Hamiltonian matrix
tial equation in Eq. (85). Using a standard symbolic solver
and solved using a phase-space Lie-algebra decoupling.
such as Mathematica, we find the solution
For more information about phase-space methods for
      continuous-variable quantum systems, see Ref. [54]. We
i −1 1 outline this decoupling methods in Sec. III B. As noted
ξ+ (t) = 2 tan tan − it − 1 , (89)
2λ− 2 above, the mapping between the phase-space solution and
the Hilbert-space solution is usually nontrivial. Some-
where we have defined times, however, the phase-space solution may yield differ-
ential equations that are simpler to solve, compared with
1 the Hilbert-space solution. It is generally difficult to say in
 2 = λ+ λ− − . (90)
4 advance which solution is the easiest to work with.

010201-14
SOLVING QUANTUM DYNAMICS. . . PRX QUANTUM 6, 010201 (2025)

C. Most general single-mode Gaussian Hamiltonian any new operators by examining the new commutator
Once we have identified two closed algebras, we can relations:
combine them to obtain more general solutions, provided 1 1
that the full algebra remains finite. For example, we can [K̂0 , â] = − â, [K̂+ , â] = − ↠,
combine the solutions for the linear and quadratic driving 2 2
(99)
terms that we derived in Secs. IV A and IV B. Doing so 1 1
[K̂0 , ↠] = ↠, [K̂− , ↠] = â.
provides us with a solution for the most general dynam- 2 2
ics of a harmonic oscillator with Gaussian terms. If the
coefficients are constant in time, we can usually compute That is, commuting the linear operators with the quadratic
them directly, as, e.g., in Ref. [66]. In general, though, the operators leaves the linear algebra invariant.
coefficients may be time dependent. We could now proceed to solve the system for the full
To solve the dynamics of a Hamiltonian with both linear algebra in Eq. (99). This means that we would have to
and quadratic interaction terms, we could start by writing compute all five multiplications by congruence according
down the full algebra and following the decoupling recipe to Eq. (51) and then solve the six resulting differential
in Sec. III. The full algebra has five unique elements, so equations simultaneously. While this is certainly possible
this necessarily involves multiplying out each of the terms (and might in this case not be too challenging), there is,
in the ansatz and solving a set of five simultaneous dif- as mentioned at the beginning of this section, an easier
ferential equation. However, we can instead make use of alternative, which involves defining a rotating frame due
the fact that we already know the solution for one of the to one subalgebra. This is similar to the notion of moving
subalgebras (e.g., the quadratic one in Sec. IV B). By then to an interaction picture and solving the dynamics there,
considering an interaction picture that rotates with one of and then returning to the laboratory frame. Commonly, the
the subalgebras, we are able to determine the effects of free-evolution Hamiltonian Ĥ 0 of the system is chosen for
the solutions on the second subalgebra. In this way, it is this but there is nothing that prevents us from choosing a
possible to partition the dynamics and solve the individ- more complex Hamiltonian, even a time-dependent one, as
ual contributions separately. It should be noted that there the rotating frame (see the Appendix A).
is no right or unique way to perform this partition and that There is no set recipe for how to best partition the
sometimes one partition works better than the other. dynamics but often one composition is easier to treat than
We start by combining the Hamiltonians in Eq. (53) and the other. In our case, we choose to focus on the linear
Eq. (77) into a single Hamiltonian with both linear and terms and how they evolve under the quadratic subalge-
quadratic terms: bra. The reason for this choice is that the quadratic algebra
leaves the linear algebra invariant, while the action of the
linear algebra on the quadratic algebra reintroduces linear
Ĥ (t) = Ĥ 0 + Ĥ L (t) + Ĥ Q (t), (96)
components into the Lie-algebra elements. We consider a
frame that rotates with the quadratic Hamiltonian Ĥ Q (t).
where we have used the subscripts L and Q to denote the The time evolution generated by the quadratic part of the
linear and quadratic terms, respectively. The Hamiltonian Hamiltonian is given by
contributions in Eq. (96) are given by
  t 
i
Ĥ 0 = ω↠â, ÛQ (t) = T exp − dt Ĥ 0 + Ĥ Q (t ) . (100)
 0
Ĥ L (t) = ωg+ (t)↠+ ωg− (t)â, (97)
Note that we have included the free-evolution Ĥ 0 in
Ĥ Q (t) = ωλ+ (t)â + ωλ− (t)â .
†2 2
Eq. (100) to complete the quadratic algebra. We already
know the solution to ÛQ (t), with the ansatz shown in
As before, ω is the angular frequency of the free mode, the Eq. (80), and the differential equations for the coefficients
g± (t) are the coefficients of the linear terms, and the λ± (t) listed in Eq. (88).
are the coefficients of the quadratic terms. Next, we consider how the Hamiltonian Ĥ L (t) evolves
The full algebra generated by the Hamiltonian in in the frame rotating with ÛQ (t). Applying this solution to
Eq. (96) is now given by ÛQ (t) in Eq. (80) to Ĥ L (t), we find that it evolves as

â, ↠, 1,  
ÛQ (t)Ĥ L (t)ÛQ (t) = Û−1
†
Q (t) g+ (t)â + g− (t)â ÛQ (t)
†
(98)
K̂0 , K̂+ , K̂− ,
= g+ (t)eiξ− K̂− eiξ0 K̂0 eiξ+ K̂+ ↠e−iξ+ K̂+ e−iξ0 K̂0 e−iξ− K̂−
with K̂0 and K̂± defined in Eq. (78). We can check that + g− (t)eiξ− K̂− eiξ0 K̂0 eiξ+ K̂+ âe−iξ+ K̂+ e−iξ0 K̂0 e−iξ− K̂− .
combining the algebras in this way does not generate (101)

010201-15
SOFIA QVARFORT and IGOR PIKOVSKI PRX QUANTUM 6, 010201 (2025)

To proceed, we first need to compute the following expres- linear algebra are given by
sions:
â, ↠, 1. (109)
−iξ0 K̂0 − 12 iξ0
e iξ0 K̂0
âe =e â,
Note that here we have not included the free-evolution
iξ0 K̂0 † −iξ0 K̂0
e â e =e
1 iξ †
2 0 â , ↠â, since that has already been included in the quadratic
(102) algebra. Just as in Sec. IV A, we make the ansatz
eiξ+ K̂+ âe−iξ+ K̂+ = â − iξ+ ↠,
†
ÛL (t) = e−iFI 1 e−iF+ â e−iF− â , (110)
eiξ− K̂− ↠e−iξ− K̂− = ↠+ iξ− â.
where the F coefficients are now different from that in
Quadratic transformations of the kind shown in the second Eq. (64), because the Hamiltonian is given by that in
two equations of Eq. (102) are also known as Bogoliubov Eq. (101). The solutions to F± in Eq. (110) are given by
transformations. They map ladder operators to linear mix-
 t  t
tures of themselves or to linear combinations of different
modes. Using the expressions in Eq. (102), we find F+ = dt ν(t ), F− = dt μ(t ) (111)
0 0

eiξ− K̂− eiξ0 K̂0 eiξ+ K̂+ âe−iξ+ K̂+ e−iξ0 K̂0 e−iξ− K̂− and
1 1 1
 t  t
= e− 2 iξ0 + ξ+ ξ− e 2 iξ0 â − iξ+ e 2 iξ0 ↠(103) FI = −i dt μ(t ) dt ν(t ). (112)
0 0

and These integrals might be challenging to perform analyt-

ically, depending on the explicit time dependence of the
eiξ− K̂− eiξ0 K̂0 eiξ+ K̂+ ↠e−iξ+ K̂+ e−iξ0 K̂0 e−iξ− K̂− coefficients, which in turn may depend on the previously
1   obtained solutions. The above expressions do, however,
= e 2 iξ0 ↠+ iξ− â . (104) capture the dynamics exactly.
We are then finally in a position to write down the full
Equation (101) therefore becomes solution to the dynamics generated by the Hamiltonian in
Eq. (96). It is given by
† 1 1
ÛQ (t)Ĥ L (t)ÛQ (t) = g+ (t)e 2 iξ0 − ig− (t)ξ+ e 2 iξ0 ↠†
Û(t) = e−iξ+ K̂+ e−iξ0 K̂0 e−iξ− K̂− e−iFI 1 e−iF+ â e−iF− â . (113)
1 1 1
+ ig+ (t)ξ− e 2 iξ0 + g− (t) e− 2 iξ0 + ξ+ ξ− e 2 iξ0 â.
By again contracting the last two exponentials as a dis-
(105) placement operator, Eq. (114) can also be written as

For convenience, we define the following coefficients Û(t) = e−iξ+ K̂+ e−iξ0 K̂0 e−iξ− K̂− e−iϕ D̂(−iF+ ), (114)
1 1
μ(t) = g+ (t)e 2 iξ0 − ig− (t)ξ+ e 2 iξ0 , where ϕ is defined in Eq. (70). This procedure, in which
1 1 1
two subalgebras are combined into a single closed algebra,
ν(t) = ig+ (t)ξ− e 2 iξ0 + g− (t) e− 2 iξ0 + ξ+ ξ− e 2 iξ0 , can in principle be repeated for more than one mode, as
(106) long as the complete algebra remains finite in terms of its
dimensions.
which allows us to write
V. NONLINEAR QUANTUM DYNAMICS: CAVITY
  OPTOMECHANICS
Û−1
Q (t) g+ (t)â + g− (t)â ÛQ (t) = μ(t)â + ν(t)â.
† †

(107) In continuous-variable systems, we may study the

dynamics of a system by deriving the equations of motion
We now state the ansatz for the linear algebra: of quantum operators such as the position or momen-
tum quadrature operators. When the Hamiltonian contains
† terms with at most quadratic quadrature operators, the
ÛL (t) = e−iFI 1 e−iF+ â e−iF− â . (108) resulting equations of motion are linear and map Gaus-
sian states to Gaussian states. If the Hamiltonian contains
We can find the solution for ÛL (t) by using the same recipe terms with products of more than two operators, the equa-
as in Sec. IV A. We know that the algebra operators of the tions of motion instead become nonlinear in terms of the

010201-16
SOLVING QUANTUM DYNAMICS. . . PRX QUANTUM 6, 010201 (2025)

operators. Compared with linear or quadratic Hamiltoni- More broadly, commuting observables are closely linked
ans, this kind of dynamics can map Gaussian states into to the notion of globally conserved physical quantities but
non-Gaussian states. There are not many nonlinear Hamil- noncommuting observables can also be studied from the
tonians that generate a Lie algebra with a finite number of perspective of Lie algebras [70].
elements. In fact, most nonlinear Hamiltonians generate an The nonlinear optomechanical Hamiltonian in Eq. (115)
infinite number of terms and must therefore be treated with is one of the few examples of a nonquadratic Hamiltonian
perturbation theory. However, there are a few examples of with a closed algebra, which allows for the dynamics to be
nonlinear dynamics that allow for analytic solutions. We solved exactly. Other examples of solutions for nonlinear
detail one of those examples here. dynamics include optomechanical systems with multiple
In cavity optomechanics, a mechanical oscillator is cou- interaction [26], as well as harmonic oscillators interacting
pled through radiation pressure to an optical mode [67]. through a cross-Kerr coupling [27].
The first-order Hamiltonian-interaction term contains a Let us apply the Lie-algebra decoupling theorem to the
product of three operators, which means that the dynamics Hamiltonian in Eq. (115). We find the following algebra
generate non-Gaussian states. The cavity optomechanical terms (this time in a Hermitian basis):
Hamiltonian reads
↠â, b̂† b̂, (↠â)2 ,
Ĥ OMS = ωc â â + ωm b̂ b̂ − g(t)â â(b̂ + b̂), (115)
† † † † (118)
↠â(b̂† + b̂), i↠â(b̂† − b̂).
where ωc is the optical oscillation frequency of the opti-
cal mode with annihilation and creation operators â, ↠, ωm After rescaling all Hamiltonian terms and time by ωm , such
is the mechanical oscillation frequency of the mechanical that ωm t → t and g(t)/ωm → g(t), we find that the ansatz
mode with annihilation and creation operators b̂ and b̂† , for the time-evolution operator Û(t) then becomes
and g(t) is the coupling strength between the optical and † † â)2 −iF ↠â(b̂† +b̂) −iF ↠â[i(b̂† −b̂)]
mechanical modes. The dynamics of this Hamiltonian with Û(t) = e−ib̂ b̂t e−iFa (â e+ −e ,
additional linear and quadratic mechanical terms have been (119)
previously solved in full generality [19–21].
We note that this Hamiltonian is the same as the linearly where we have transformed into a frame  that rotates with
driven quantum harmonic oscillator explored in Sec. IV A, the free optical evolution exp −i↠âtωc .
except that the driving term is now multiplied with the The F coefficients in Eq. (119) are given by the follow-
operator ↠â. If we proceed to map out the algebra of this ing integrals [19,20]:
system by taking the commutator between the free evolu-  
tion of the mechanical mode and the interaction term in t   t  
Eq. (115), we find the following term: Fa = 2 dt g(t ) sin t dt g(t ) cos t ,
0 0
 t  
[b̂† b̂, ↠â(b̂† + b̂)] = ↠â(b̂† − b̂). (116) F+ = − dt g(t ) cos t , (120)
0

Then, commuting this new term with the original interac- t  
tion term, we find F− = dt g(t ) sin t .
0

[â â(b̂ + b̂), â â(b̂ − b̂)] = 2(â â) .

† † † † † 2
(117) For a constant optomechanical coupling g(t) ≡ g0 , the
integrals in Eq. (120) evaluate to
The operator (↠â)2 corresponds to a new algebra element.
It is a self-Kerr interaction that can also be identified by 1 2
diagonalizing the Hamiltonian using a polaron transform Fa = g [sin(2t) − 2t] ,
2 0 (121)
[68,69]. In this context, it is also known as a Casimir F+ = − g0 sin(t), F− =g0 [cos(t) − 1] ,
element. A Casimir element is an element that is part of
the algebra but that commutes with all other operators which is equivalent to the previously obtained solutions
in the algebra. A trivial example is the identity opera- [68,69] up to the ordering of the terms in Eq. (119).
tor 1, which arises as part of the linear algebra shown in
Eq. (54) from which we take the following commutator:
VI. DECOUPLING THEOREM FOR
[↠+ â, i(↠− â)] = 2i1. Since the identity operator com-
OPEN-SYSTEM DYNAMICS
mutes with every other element in the algebra, it can be
considered a Casimir element. In the case of the self-Kerr In all previous sections, we have considered quantum
nonlinearity (↠â)2 , however, we cannot discard it, since systems isolated from their environment. However, in all
it usually leads to observable changes in a quantum state. realistic quantum systems, there is at least some coupling

010201-17
SOFIA QVARFORT and IGOR PIKOVSKI PRX QUANTUM 6, 010201 (2025)

to the surroundings. While the Lie-algebra decoupling To vectorize the Lindblad equation, we need to make use
method might at first seem useful only for closed dynam- of the following relation (for the derivation, see Ref. [74])
ics, it can be readily applied toward solving open-system
dynamics as well. |ABC = (Â ⊗ ĈT )|B. (125)
There are many ways to model open-system dynamics.
Here, we focus on the Gorini-Kossakowski-Sudarshan- Our goal is to later replace B̂ by the density matrix in the
Lindblad equation (referred to from here on as the Lind- Lindblad equation, such that it can be written as a linear
blad equation), which is the most general Markovian matrix equation.
master equation [35,36], Finally, we note that the expectation value for the gen-
eral operator  and the state ˆ is given in the vectorized
2 −1   language as
N
˙ˆ = −i[Ĥ , ]
ˆ +
1
hnm L̂n ˆ L̂†m − {L̂†m L̂n , }
ˆ , (122)  
2  = Tr ˆ = A† |. (126)
n,m=1
This will allow us to compute various quantities of interest
once we have solved the dynamics.
where ˆ is the density matrix of a quantum state, Ĥ is the
Hamiltonian operator, L̂n is a Lindblad operator, and {·, ·} B. Vectorizing the Lindblad equation
denotes the anticommutator.
In order to apply the Lie-algebra decoupling method As we have seen, vectorization allows us to turn matri-
to the Lindblad equation in Eq. (122), we must first ces into vectors. Crucially, it also allows us to transform
rewrite it so that we can state the ansatz with products of superoperators into matrices. In this work, we denote the
exponentials. To do so, we vectorize the state. vectorized density matrix ˆ by | and the free-state
evolution is subsequently written as

A. Introduction to vectorization Û(t)ˆ 0 Û† (t) → Û(t) ⊗ Û∗ (t)|0 . (127)
Intuitively, the vectorization procedure turns an opera- Note that here we take the complex conjugate rather than
tor or a density matrix into a vector by stacking its rows the full conjugate transpose of Û(t), as mandated by the
and columns. The choice of row or column depends on vectorization mapping we have chosen. The tensor prod-
convention and both are equally valid. Here, we intro- uct is used to differentiate between the left-hand and
duce the vectorization procedure for linear operators that right-hand multiplication of Û(t).
act on the Hilbert space and show how the vectorized To apply the vectorization to the Lindblad equation, we
Lindblad equation is derived. We also refer the reader use the identity in Eq. (125) on all terms of the Lindblad
to the introduction to vectorization in Refs. [37,71,72]. equation in Eq. (122). Each term and its corresponding
See also Ref. [73] and in the Supplemental Material of vectorized version can be found in Table I. Where only two
Ref. [74], the notation of which we follow closely. Alterna- operators were multiplied, we have inserted the identity
tives to vectorization include rewriting the density matrix operator to obtain a product of three operators. As a result,
in Louville-Bloch form [47]. Eq. (122) can be written in the vectorized language as
We start by considering a generic operator  that acts on
the Hilbert space H. Given an orthonormal basis {|i} in d
| = L̂(t)|. (128)
H, the operator  can be written as dt

 Here, we write L̂(t) as:

 = i|  |j  |ij | . (123)
ij L̂ = L̂H + L̂L , (129)

TABLE I. Vectorized analogues of terms in the Lindblad

We then assign a vector to this operator by flipping one of
equation [Eq. (122)].
the bras into a ket, such that
Operator product Vectorized analogue

|A = i| Â |j  |i |j  . (124) ˆ Ĥ (1 ⊗ Ĥ T )|
ij Ĥ ˆ (Ĥ ⊗ 1)|
† †T
L̂n ˆ L̂m (L̂n ⊗ L̂m )|
† †
L̂m L̂n ˆ (L̂m L̂n ⊗ 1)|
That is, we define |A by stacking every row in the matrix † †
ˆ L̂m L̂n (1 ⊗ (L̂m L̂n )T )|
Â, which defined through its elements Aij = i| Â |j .

010201-18
SOLVING QUANTUM DYNAMICS. . . PRX QUANTUM 6, 010201 (2025)

where (according to the terms listed in Table I) L̂H is algebra. Therefore, to the best of our knowledge, there are
the unitary (Hamiltonian) contribution given by L̂H = not many examples of analytically solvable open-system
dynamics. Here, we outline the case of dissipation and
−i Ĥ (t) ⊗ 1 − 1 ⊗ Ĥ T (t) and L̂L contains the nonuni- thermalization, which are commonly used to describe a
tary part variety of quantum systems. For certain noise baths, a
2 −1 quantum system will inevitably thermalize with its sur-
N
hnm rounding temperature. Such a process involves both a
L̂L = 2L̂n ⊗ L̂†T
m − L̂m L̂n ⊗ 1 + 1 ⊗ (L̂m L̂n )
† † T
.
dissipative process and a heating process, which ensures
n,m=1
2
that the equilibrium temperature is achieved.
(130) We can model the interplay of heating and dissipa-
tion by assuming that the Lindblad equation contains two
We note the appearance of transposed operators in √
Eq. (130), which follow from our choice of the vectoriza- different Lindblad operators, namely, L̂1 = κnth ↠and
√
tion mapping. We may simplify the expression by adopting L̂2 = κ(1 + nth )â. The resulting Lindblad equation reads
a real basis, such as the Fock basis, where L̂ and L̂† have (again setting  = 1)
exclusively real entries. This means that the transposition
operation is equivalent to taking the Hermitian conjugate, d κnth  † 
† (t)
ˆ = −i[Ĥ , (t)]
ˆ + 2â (t)â − {â↠, (t)}
ˆ
which, e.g., allows us to write L̂Ti = L̂i . This will greatly dt 2
simplify our calculations but we might have to be careful κ(1 + nth )  
when we compute quantities defined in a complex basis, + 2â(t)â
ˆ †
− {↠â, (t)} . (134)
2
such as coherent states. Considering coherent states in the
Fock basis, however, is fine. To solve this equation, we apply the Lie-algebra decou-
The formal solution to the Lindblad equation in pling theorem. The derivation has also been carried out
Eq. (128) in the vectorized language reads in Ref. [75] (where vectorization is referred to as the
thermofield method). The vectorized Louvillian that cor-
|(t) = Ŝ (t)|0 , (131) responds to the above Lindblad equation (where we have
where |0  is the vectorized form of the initial state ˆ 0 and also removed the transposes by assuming that we are
working in a real basis) is
Ŝ (t) is the time-ordered exponential of L̂(t):
 t 
←
− L̂ = −i Ĥ ⊗ 1 − 1 ⊗ Ĥ
Ŝ (t) = T exp dt L̂(t ) . (132)
0
κ(nth + 1)  
This is a key expression that captures both the unitary and + 2â ⊗ â − ↠â ⊗ 1 − 1 ⊗ ↠â
2
the nonunitary evolution. κnth  † 
Now that we have written the dynamical operator in + 2â ⊗ ↠− â↠⊗1 − 1 ⊗ â↠. (135)
2
exponential form, we can make the following ansatz:
Based on this expression, we can identify the following Lie
Ŝ (t) = Ŝ1 (t)Ŝ2 (t) . . . ŜN , (133)
algebra that generates the evolution of the system
where each term is defined as Ŝj = eDj L̂j . Here, Dj is a
dimensionless time-dependent coefficient and L̂j is a Lie- K̂+ =↠⊗ ↠, K̂− =â ⊗ â,
algebra element. 1 † 
Crucially, the Lie-algebra decoupling theorem only K̂3 = â â ⊗ 1 + 1 ⊗ ↠â + 1 , (136)
2
requires us to take the inverse of operators. For the expo-
K̂0 = ↠â ⊗ 1 − 1 ⊗ ↠â.
nential map eA , its inverse is always defined as e−A , which
allows us to follow the steps outlined in Sec. III C and
derive differential equations for the Dj coefficients. This We may immediately note that this algebra is the two-mode
allows us to consider and solve open dynamics where the version of the quadratic algebra identified in Sec. IV B.
Hamiltonian and Lindbladian generate a finite-dimensional By then assuming that the Hamiltonian is just a free har-
Lie algebra. We demonstrate such a case in Sec. VII. monic oscillator Ĥ = ω0 ↠â, we can write the Louvillian
in terms of the algebra elements
VII. OPEN DYNAMICS: A THERMALIZING
HARMONIC OSCILLATOR L̂ = −iω0 K̂0 + κ(nth + 1)K̂− +κnth K̂+
Linear Lindblad operators give rise to quadratic terms in κ
− κ(2nth + 1)K̂3 + . (137)
the Lindblad equation, which means that we have finite Lie 2

010201-19
SOFIA QVARFORT and IGOR PIKOVSKI PRX QUANTUM 6, 010201 (2025)

The ansatz for the evolution of the system then reads First, if the algebra of the Hamiltonian does not close,
it becomes impossible to derive analytic expressions for
Ŝ (t) = eγ0 eD0 K̂0 eD+ K̂+ eD3 K̂3 eD− K̂− , (138) the associated differential equations. In these cases, one
instead has an infinite product of exponentials like the
where D0 = −iωt and where the symbols D± and D0 give Zassenhaus formula, or the Trotter expansion. However,
rise to nonunitary dynamics, compared with the coeffi- it may still be feasible to leverage the method pertur-
cients that we have found for the quadratic Hamiltonian batively, extending it to higher orders of the interaction
in Sec. IV B. terms. An example of that is the Jaynes-Cummings model
By then applying the decoupling theorem, we can derive [76]. Second, there are situations in which the method
solutions for the coefficients that are very similar to those can technically be applied but may offer limited practi-
identified in Eq. (95). By defining cal advantage. A notable example is the still-Gaussian
case with an interaction term of the form (↠+ â)(b̂† + b̂).
γ+ =κnth t, γ− =κ(nth + 1)t, Commuting this term with the free evolution of either field
1 (139) generates ten distinct terms. While it is straightforward
γ3 = −(γ+ +γ− ), γ0 = (γ− −γ+ ), although tedious to formulate the corresponding differ-
2
ential equations, depending on the time-dependent coef-
we find that the coefficients in Eq. (138) are given by ficients, solving them can be highly challenging. Lastly,
another demanding scenario arises in many-body inter-
D0 = −iωt, actions involving multiple modes. In such cases, the Lie
2γ± sinh(φ) algebra of a multimode system can quickly become too
D± = , complex to manage. Nevertheless, certain patterns or sym-
(2φ cosh(φ) − γ3 sinh(φ)) (140) metries might emerge, enabling the application of the
1 method in specific contexts. Despite these challenges, par-
D3 = 2
,
γ3 titioning the dynamics into solvable and unsolvable com-
cosh(φ) − 2φ sinh(φ)
ponents, where the solvable components are treated with
the Lie-algebra method, can still yield valuable insights,
where we also have defined particularly when combined with techniques such as time-
dependent perturbation theory.
γ32
φ2 = − γ+ γ− . (141) To conclude, the goal of this tutorial has been to pro-
4 vide valuable insights on the dynamics of quantum systems
The resulting dynamics give rise to a state that gradually and practical guidance on how to solve the time evolu-
either dissipates or thermalizes toward a thermal state with tion of open and closed quantum systems. The Lie-algebra
occupation number nth . decoupling method offers a powerful and versatile tool set
for tackling complex quantum systems with driving terms
VIII. SUMMARY of arbitrary time dependence. The methods presented here
apply to a number of Hamiltonians and physical systems,
In this tutorial, we have provided a pedagogical intro- thereby paving the way for approaching many intriguing
duction to solving the dynamics of quantum systems future research avenues in quantum physics.
using a Lie-algebra decoupling method. As a demon-
stration of the method, we have considered both uni-
tary and nonunitary dynamics. The examples we have
covered have included a quantum harmonic oscillator
ACKNOWLEDGMENTS
with a linear and quadratic time-dependent single-mode
Hamiltonian interaction term and a nonlinear Hamiltonian We thank Yuefei Liu, Suocheng Zhao, Sreenath K.
found in the context of cavity optomechanics, as well as Manikandan, and David Edward Bruschi for fruitful dis-
one of the most common noise models, which include cussions and comments. S.Q. is funded in part by the
both dissipation and thermalization. For each example, Wallenberg Initiative on Networks and Quantum Informa-
we have applied the method step by step and derived tion (WINQ) and in part by the Marie Skłodowska-Curie
the exact differential equations that govern the dynam- Action Individual Fellowship program Nonlinear Optome-
ics. The models that we have considered arise in many chanics for Verification, Utility, and Sensing (NOVUS)
situations in quantum optics and related fields, with appli- under Grant No. 101027183. Nordita is funded in part
cations for both quantum technologies and fundamental by NordForsk. I.P. acknowledges support by the Swedish
physics. Research Council under Grant No. 2019-05615, the Euro-
As with all methods, there are scenarios in which the pean Research Council under Grant No. 742104 and The
Lie-algebra decoupling approach proves less effective. Branco Weiss Fellowship—Society in Science.

010201-20
SOLVING QUANTUM DYNAMICS. . . PRX QUANTUM 6, 010201 (2025)

APPENDIX: PARTITIONING QUANTUM Schrödinger picture shown in Eq. (A2), we find

DYNAMICS
d d †
In this appendix, we provide a derivation of the time i |ψ(t)I = Û0 (t) |ψ(t)S
partitioning of the evolution operators, which we have dt dt
made use of in Sec. IV C. We start by considering unitary i †
− Û0 (t) Ĥ 0,S (t) + Ĥ 1,S (t) |ψ(t)S .
dynamics but then we generalize the results to nonunitary 
dynamics in the vectorized language. (A6)

†
1. Unitary dynamics In addition, from the definition of Û0 (t) in Eq. (A4), its
derivative is given by
This proof follows the usual derivation of the interac-
tion picture; however, we show that any partition of the  †  †
d † d i
Hamiltonian will do, even one where the two parts are time Û (t) = Û0 (t) = − Ĥ 0,S (t)Û0 (t)
dependent. dt 0 dt 
We start by writing down the Hamiltonian in the i †
Schrödinger picture, where the S indices denote a quantity = Û0 (t)Ĥ 0,S (t), (A7)

in the Schrödinger picture. In contrast, the I indices denote
the interaction picture. Let us start with the Hamiltonian which follows from Leibniz’s integral rule and the fact
that the adjoint operation commutes with the differentiation
Ĥ S (t) = Ĥ 0,S (t) + Ĥ 1,S (t). (A1) operator. Inserting Eq. (A7) into Eq. (A5), we find that

d i †
In the standard derivation, Ĥ 0,S (t) is often taken to the i |ψ(t)I = Û0 (t)Ĥ 0,S (t) |ψ(t)S
dt 
free (time-independent) evolution of the system. While this i †
is often a convenient choice (especially when performing − Û0 (t) Ĥ 0,S (t) + Ĥ 1,S (t) |ψS

approximations, such as the rotating-wave approximation),
any convenient partitioning is allowed. ≡ Ĥ 1,I (t) |ψ(t)I , (A8)
We proceed by writing down the Schrödinger equation
in the Schrödinger picture: where we have defined the Hamiltonian in the interaction
picture as
d
i |ψ(t)S = Ĥ 0,S (t) + Ĥ 1,S (t) |ψ(t)S . (A2) †
Ĥ 1,I (t) = Û0 (t)Ĥ 1,S (t)Û0 (t). (A9)
dt

Next, we define the state in the interaction picture, The formal solution to Eq. (A8) is, as usual,

† |ψ(t)I = Û1 (t) |ψ(0)I , (A10)

|ψ(t)I = Û0 (t) |ψ(t)S , (A3)
where we have defined Û1 (t) as
where Û0 (t) is given by the time-ordered exponential over
  t 
Ĥ 0,S i †
Û1 (t) = T exp − dt Û0 (t)Ĥ 1,S (t)Û0 (t) . (A11)

   0
i t
Û0 (t) = T exp − dt Ĥ 0,S (t ) . (A4) To complete our derivation, we ask what the full state in
 0
the Schrödinger picture looks like. Multiplying Eq. (A3)
by Û0 (t) on both sides, we find
We then consider the Schrödinger equation in the interac-
tion picture. We find
|ψ(t)S = Û0 (t) |ψ(t)I = Û0 (t)Û1 (t) |ψ(0) . (A12)

d d † † d We have shown that any arbitrary partition of the Hamilto-

i |ψ(t)I = Û0 (t) |ψ(t)S + Û0 (t) |ψ(t)S . (A5)
dt dt dt nian, and thereby the time evolution, allows us to pairwise
multiply the state by the time-evolution operators. There
We note the appearance of the derivative of |ψ(t)S . Insert- is no unique partition and it is often a question of what is
ing the right-hand side of the Schrödinger equation in the most convenient.

010201-21
SOFIA QVARFORT and IGOR PIKOVSKI PRX QUANTUM 6, 010201 (2025)

2. Nonunitary dynamics As a result, the derivative of Ŝ0−1 (t) is

The question now arises of whether we can perform the
same partition for nonunitary evolution. Such a partition d −1
Ŝ (t) = −Ŝ0−1 (t)L̂0,S (t). (A21)
could be useful when both unitary and nonunitary terms dt 0
contribute to the dynamics. To prove that a partition is pos-
sible, we make use of the vectorization notation introduced Putting everything together, we insert the result in
in Sec. VI. Let the full Lindbladian (which includes both Eq. (A21) into Eq. (A17). We also note that the first
unitary and nonunitary terms) be term on the right-hand side of Eq. (A17) is given by the
Lindblad equation. Inserting this, we find
L̂(t) = L̂0 (t) + L̂1 (t), (A13)
d
|(t)I = −Ŝ0−1 (t)L̂0,S |(t)S
where L̂0 (t) is the term that we wish to solve first, and then dt
apply to L̂1 (t). Then, the vectorized Lindblad equation in + Ŝ0−1 (t) L̂0,S + L̂1,S |(t)S . (A22)
the Schrödinger picture, which we denote by the subscript
S (see the main text for the derivation), reads
Canceling the terms and defining
d
|(t)S = L̂0,S (t) + L̂1,S (t) |(t)S . (A14)
dt L̂1,I (t) = e−L̂0,S t L̂1,S eL̂0,S t , (A23)
We then define the vectorized state of the system in the we have
interaction picture as
d
|(t)I = Ŝ0 (t)−1 |(t)S , (A15) |(t)I = L̂1,I (t)|(t)I . (A24)
dt

where Ŝ0 (t) is defined as The formal solution to Eq. (A24) is

 t 
Ŝ0 (t) = T exp dt L̂0 (t ) . (A16) |(t)I = Ŝ1 (t)|(0)I , (A25)
0

To determine the evolution of |(t)I in the interaction where Ŝ1 is defined as

picture, we consider the Lindblad equation, shown in  
←
− t
Eq. (128). Taking the derivative of |(t)I , we find Ŝ1 (t) = T exp dt e−L̂0,S t L̂1,S (t )eL̂0,S t . (A26)
0
d d
|(t)I = Ŝ −1 (t)|(t)S Finally, we wish to return to the Schrödinger picture.
dt dt 0
d d Multiplying (A15) on both sides by Ŝ0 (t), we find
= Ŝ0−1 (t)|(t)S + Ŝ0−1 (t) |(t)S . (A17)
dt dt
|(t)S = Ŝ0 (t)|(t)I = Ŝ0 (t)Ŝ1 (t)|(0). (A27)
To evaluate the first term on the right-hand side of
Eq. (A17), we have to compute the derivative of an This proves that nonunitary time evolution can also be
inverted matrix. Using the identity 1 = Ŝ0−1 (t)Ŝ0 (t), which partitioned into two arbitrary terms, just as is the case
is well defined for the exponential map, we find for unitary dynamics. Note that we have not made any
assumptions about the Hermiticity of either Ŝ0 (t) or Ŝ1 (t),
d d d
1 = 0 = Ŝ0−1 (t)Ŝ0 (t) + Ŝ0−1 (t) Ŝ0 (t). (A18) but that everything follows from the fact that the inverse
dt dt dt is well defined for exponential maps. One useful strategy
Rearranging Eq. (A18), we find could be, e.g., to partition the dynamics into unitary and
nonunitary dynamics, such that Ŝ0 (t) contains the unitary
d −1 d part and ŜI (t) contains the nonunitary part rotating in the
Ŝ0 (t) = −Ŝ0−1 (t) Ŝ0 (t)Ŝ0−1 (t). (A19)
dt dt frame defined by Ŝ0 (t).
The derivative of Ŝ0 (t) again follows from Leibniz’s inte-
gral rule, which gives us
[1] D. Dong and I. R. Petersen, Quantum control theory and
d
Ŝ0 (t) = L̂0,S (t)Ŝ0 (t). (A20) applications: A survey, IET Control Theory Appl. 4, 2651
dt (2010).

010201-22
SOLVING QUANTUM DYNAMICS. . . PRX QUANTUM 6, 010201 (2025)

[2] A. Blais, J. Gambetta, A. Wallraff, D. I. Schuster, S. M. optomechanical systems: Interplay of mechanical squeez-
Girvin, M. H. Devoret, and R. J. Schoelkopf, Quantum- ing and non-Gaussianity, J. Phys. A: Math. Theor. 53,
information processing with circuit quantum electrodynam- 075304 (2020).
ics, Phys. Rev. A 75, 032329 (2007). [22] F. Schneiter, S. Qvarfort, A. Serafini, A. Xuereb, D. Braun,
[3] C. L. Degen, F. Reinhard, and P. Cappellaro, Quantum D. Rätzel, and D. E. Bruschi, Optimal estimation with
sensing, Rev. Mod. Phys. 89, 035002 (2017). quantum optomechanical systems in the nonlinear regime,
[4] M. Maggiore, A generalized uncertainty principle in quan- Phys. Rev. A 101, 033834 (2020).
tum gravity, Phys. Lett. B 304, 65 (1993). [23] S. Qvarfort, A. D. K. Plato, D. E. Bruschi, F. Schneiter,
[5] A. Kempf, G. Mangano, and R. B. Mann, Hilbert space rep- D. Braun, A. Serafini, and D. Rätzel, Optimal estima-
resentation of the minimal length uncertainty relation, Phys. tion of time-dependent gravitational fields with quan-
Rev. D 52, 1108 (1995). tum optomechanical systems, Phys. Rev. Res. 3, 013159
[6] I. Pikovski, M. R. Vanner, M. Aspelmeyer, M. Kim, and (2021).
Č. Brukner, Probing Planck-scale physics with quantum [24] S. Qvarfort, M. R. Vanner, P. F. Barker, and D. E. Bruschi,
optics, Nat. Phys. 8, 393 (2012). Master-equation treatment of nonlinear optomechanical
[7] F. Armata, L. Latmiral, I. Pikovski, M. R. Vanner, Č. systems with optical loss, Phys. Rev. A 104, 013501 (2021).
Brukner, and M. Kim, Quantum and classical phases in [25] Y. Ling, S. Qvarfort, and F. Mintert, Fast optomechanical
optomechanics, Phys. Rev. A 93, 063862 (2016). photon blockade, Phys. Rev. Res. 5, 023148 (2023).
[8] L. Latmiral, F. Armata, M. G. Genoni, I. Pikovski, and [26] D. E. Bruschi, Time evolution of coupled multimode and
M. S. Kim, Probing anharmonicity of a quantum oscilla- multiresonator optomechanical systems, J. Math. Phys. 60,
tor in an optomechanical cavity, Phys. Rev. A 93, 052306 062105 (2019).
(2016). [27] D. E. Bruschi, Time evolution of two harmonic oscilla-
[9] H. Zassenhaus, Ein Verfahren, jeder endlichenp-Gruppe tors with cross-Kerr interactions, J. Math. Phys. 61, 032102
einen Lie-Ring mit der Charakteristikp zuzuordnen, Abh. (2020).
Math. Semin. Univ. Hambg. 13, 200 (1939). [28] S. K. Manikandan and S. Qvarfort, Optimal quantum para-
[10] W. Magnus, On the exponential solution of differential metric feedback cooling, Phys. Rev. A 107, 023516 (2023).
equations for a linear operator, Commun. Pure Appl. Math. [29] L. R. Bakker, M. S. Bahovadinov, D. V. Kurlov, V. Gritsev,
7, 649 (1954). A. K. Fedorov, and D. O. Krimer, Driven-dissipative time
[11] K. Kumar, On expanding the exponential, J. Math. Phys. 6, crystalline phases in a two-mode bosonic system with Kerr
1928 (1965). nonlinearity, Phys. Rev. Lett. 129, 250401 (2022).
[12] R. Wilcox, Exponential operators and parameter differenti- [30] J. N. Bandyopadhyay and J. Thingna, Floquet engineer-
ation in quantum physics, J. Math. Phys. 8, 962 (1967). ing of Lie algebraic quantum systems, Phys. Rev. B 105,
[13] M. Suzuki, Generalized Trotter’s formula and systematic L020301 (2022).
approximants of exponential operators and inner deriva- [31] Y. Gao, N. Yu, S. Wang, and G. Wang, Simulation of mixed
tions with applications to many-body problems, Commun. quantum Rabi model and its applications on generation of
Math. Phys. 51, 183 (1976). squeezed cat state, Int. J. Theor. Phys. 61, 1 (2022).
[14] M. Suzuki, On the convergence of exponential opera- [32] H.-C. Fu and Z.-R. Gong, Exact solution for non-
tors—the Zassenhaus formula, BCH formula and system- Markovian master equation using hyper-operator approach,
atic approximants, Commun. Math. Phys. 52, 193 (1977). Commun. Theor. Phys. 71, 1089 (2019).
[15] M. Suzuki, Relationship between d-dimensional quan- [33] L. Teuber and S. Scheel, Solving the quantum master
tal spin systems and (d + 1)-dimensional Ising Systems: equation of coupled harmonic oscillators with Lie-algebra
Equivalence, critical exponents and systematic approxi- methods, Phys. Rev. A 101, 042124 (2020).
mants of the partition function and spin correlations, Prog. [34] L. Bakker, V. Yashin, D. Kurlov, A. Fedorov, and V.
Theor. Phys. 56, 1454 (1976). Gritsev, Lie-algebraic approach to one-dimensional transla-
[16] M. Suzuki, Decomposition formulas of exponential opera- tionally invariant free-fermionic dissipative systems, Phys.
tors and Lie exponentials with some applications to quan- Rev. A 102, 052220 (2020).
tum mechanics and statistical physics, J. Math. Phys. 26, [35] G. Lindblad, On the generators of quantum dynamical
601 (1985). semigroups, Commun. Math. Phys. 48, 119
[17] N. Hatano and M. Suzuki, in Quantum Annealing and (1976).
Other Optimization Methods (Springer, Springer Berlin, [36] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Com-
Heidelberg, 2005), p. 37. pletely positive dynamical semigroups of N-level systems,
[18] J. Wei and E. Norman, Lie algebraic solution of linear J. Math. Phys. 17, 821 (1976).
differential equations, J. Math. Phys. 4, 575 (1963). [37] J. A. Gyamfi, Fundamentals of quantum mechanics in
[19] D. E. Bruschi and A. Xuereb, “Mechano-optics”: An Liouville space, Eur. J. Phys. 41, 063002 (2020).
optomechanical quantum simulator, New J. Phys. 20, [38] G. D. Mahan, Many-Particle Physics (Springer New York,
065004 (2018). NY, 2000), Chapter 4
[20] S. Qvarfort, A. Serafini, A. Xuereb, D. Rätzel, and D. [39] P. Moosavi, Exact Dirac–Bogoliubov–de Gennes dynamics
E. E. Bruschi, Enhanced continuous generation of non- for inhomogeneous quantum liquids, Phys. Rev. Lett. 131,
Gaussianity through optomechanical modulation, New J. 100401 (2023).
Phys. 21, 055004 (2019). [40] A. Klein and E. Marshalek, Boson realizations of Lie alge-
[21] S. Qvarfort, A. Serafini, A. Xuereb, D. Braun, D. bras with applications to nuclear physics, Rev. Mod. Phys.
Rätzel, and D. E. Bruschi, Time-evolution of nonlinear 63, 375 (1991).

010201-23
SOFIA QVARFORT and IGOR PIKOVSKI PRX QUANTUM 6, 010201 (2025)

[41] G. Campolieti and B. Sanctuary, The Wei-Norman [59] J. Gieseler, B. Deutsch, R. Quidant, and L. Novotny,
Lie-algebraic technique applied to field modulation in Subkelvin parametric feedback cooling of a laser-trapped
nuclear magnetic resonance, J. Chem. Phys. 91, 2108 nanoparticle, Phys. Rev. Lett. 109, 103603 (2012).
(1989). [60] P. Kinsler and P. D. Drummond, Quantum dynamics of the
[42] A. Rau and K. Unnikrishnan, Evolution operators and wave parametric oscillator, Phys. Rev. A 43, 6194 (1991).
functions in a time-dependent electric field, Phys. Lett. A [61] G. Chiribella, G. M. D’Ariano, and P. Perinotti, Applica-
222, 304 (1996). tions of the group SU(1,1) for quantum computation and
[43] B. Shadwick and W. Buell, Unitary integration: A numer- tomography, Laser Phys. 16, 1572 (2006).
ical technique preserving the structure of the quantum [62] N. Quesada, L. Helt, M. Menotti, M. Liscidini, and J. Sipe,
Liouville equation, Phys. Rev. Lett. 79, 5189 (1997). Beyond photon pairs—nonlinear quantum photonics in the
[44] A. R. P. Rau, Unitary integration of quantum Liouville- high-gain regime: A tutorial, Adv. Opt. Photonics 14, 291
Bloch equations, Phys. Rev. Lett. 81, 4785 (1998). (2022).
[45] D. B. Uskov and A. R. P. Rau, Geometric phase for N - [63] J. Guerrero and M. Berrondo, Semiclassical interpretation
level systems through unitary integration, Phys. Rev. A 74, of Wei–Norman factorization for SU (1, 1) and its related
030304 (2006). integral transforms, J. Math. Phys. 61, 082107 (2020).
[46] S. Vinjanampathy and A. R. P. Rau, Bloch sphere-like con- [64] F. Calogero, Variable Phase Approach to Potential Scatter-
struction of SU(3) Hamiltonians using unitary integration, ing (Elsevier, Sezione di Roma, Rome, Italy, 1967).
J. Phys. A: Math. Theor. 42, 425303 (2009). [65] S. M. Barnett and P. M. Radmore, Methods in Theoretical
[47] A. R. P. Rau and R. Wendell, Embedding dissipation and Quantum Optics (Oxford University Press, Oxford, 2002),
decoherence in unitary evolution schemes, Phys. Rev. Lett. Vol. 15.
89, 220405 (2002). [66] I. Pikovski, Macroscopic Quantum Systems and Gravita-
[48] W. Fulton and J. Harris, Representation Theory: A First tional Phenomena, Department of Physics, Ph.D. thesis,
Course (Springer Science+Business Media, New York, University of Vienna, 2014.
2004), Vol. 129. [67] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Cavity
[49] J. B. Gutowski, in DAMTP, Centre for Mathematical Sci- optomechanics, Rev. Mod. Phys. 86, 1391 (2014).
ences (University of Cambridge, 2007). [68] S. Bose, K. Jacobs, and P. Knight, Preparation of nonclas-
[50] A field is a fundamental algebraic structure in the form sical states in cavities with a moving mirror, Phys. Rev. A
of a set, where addition, subtraction, multiplication, and 56, 4175 (1997).
division are defined. [69] S. Mancini, V. Man’ko, and P. Tombesi, Ponderomotive
[51] S. Qvarfort, Ph.D. thesis, Quantum metrology with optome- control of quantum macroscopic coherence, Phys. Rev. A
chanical systems in the nonlinear regime, Department 55, 3042 (1997).
of Physics and Astronomy, University College London, [70] N. Yunger Halpern and S. Majidy, How to build Hamil-
[Link] 2020. tonians that transport noncommuting charges in quantum
[52] E. A. Coddington, N. Levinson, and T. Teichmann, Theory thermodynamics, npj Quantum Inf. 8, 1 (2022).
of ordinary differential equations, Phys. Today 9, 18 (1956). [71] A. Gilchrist, D. R. Terno, and C. J. Wood, Vectorization of
[53] S. L. Braunstein and P. Van Loock, Quantum informa- quantum operations and its use, ArXiv:0911.2539.
tion with continuous variables, Rev. Mod. Phys. 77, 513 [72] M. Am-Shallem, A. Levy, I. Schaefer, and R. Kosloff,
(2005). Three approaches for representing Lindblad dynamics by
[54] A. Serafini, Quantum Continuous Variables: A Primer of a matrix-vector notation, ArXiv:1510.08634.
Theoretical Methods (CRC Press, Boca Raton, 2017). [73] G. D’Ariano, P. Lo Presti, and M. Sacchi, Bell measure-
[55] D. E. Bruschi, A. R. Lee, and I. Fuentes, Time evolution ments and observables, Phys. Lett. A 272, 32 (2000).
techniques for detectors in relativistic quantum informa- [74] S. Alipour, M. Mehboudi, and A. T. Rezakhani, Quantum
tion, J. Phys. A: Math. Theor. 46, 165303 (2013). metrology in open systems: Dissipative Cramér-Rao bound,
[56] D. F. Walls and G. J. Milburn, Quantum Optics (Springer- Phys. Rev. Lett. 112, 120405 (2014).
Verlag Berlin, Heidelberg, 2008). [75] S. Chaturvedi and V. Srinivasan, Solution of the master
[57] A. Rau and D. Uskov, Effective Hamiltonians in quantum equation for an attenuated or amplified nonlinear oscilla-
physics: Resonances and geometric phase, Phys. Scr. 74, tor with an arbitrary initial condition, J. Mod. Opt. 38, 777
C31 (2006). (1991).
[58] Note that this is not a global phase term, since FI is not [76] D. Braak, Integrability of the Rabi model, Phys. Rev. Lett.
necessarily real. 107, 100401 (2011).

010201-24