Fast thermal state preparation beyond native interactions

Alexander van Lomwel,1 Paul M. Schindler,2 Modesto Orozco-Ruiz,1

Marin Bukov,2 Nguyen H. Le,1, 3 and Florian Mintert1
1
Blackett Laboratory, Imperial College London, SW7 2AZ, United Kingdom
2
Max Planck Institute for the Physics of Complex Systems, MPI-PKS, Nöthnitzerstr. 38, 01187 Dresden, Germany
3
Qolumbus SAS, 78 Boulevard de la République, 92100 Boulogne-Billancourt, France
While questions on quantum simulation of ground state physics are mostly focussed on the re-
alization of effective interactions, most work on quantum simulation of thermal physics explores
the realization of dynamics towards a thermal mixed state under native interactions. Many open
questions that could be answered with quantum simulations, however, involve thermal states with
respect to synthetic interactions.
We present a framework based solely on unitary dynamics to design quantum simulations for thermal
states with respect to Hamiltonians that include non-native interactions, suitable for both present-
day digital and analogue devices. By classical means, our method finds the control sequence to
arXiv:2601.04810v1 [quant-ph] 8 Jan 2026

reach a target thermal state for system sizes well out of reach of state-vector or density-matrix
control methods, even though quantum hardware is required to explicitly simulate the thermal state
dynamics. With the illustrative example of the cluster Ising model that includes non-native three-
body interactions, we find that required experimental resources, such as the total evolution time,
are independent of temperature and criticality.

I. INTRODUCTION of Hamiltonians with non-native interactions–such as

three-body interactions. Many of the actual physics
How to prepare finite-temperature quantum states on problems that we aim to address in quantum simulation
quantum simulation devices is a foundational question require the ability to both engineer synthetic interactions
with wide ranging implications for the applicability of and to ensure evolution to the desired thermal state in
quantum simulators. It gives controlled access to equilib- the presence of these interactions [23–25].
rium thermodynamics directly from microscopic Hamil- In this work, we present a thermal state prepara-
tonians, and indirectly to non-equilibrium thermodynam- tion framework that is compatible with synthetic inter-
ics via Jarzynski’s equality [1], provides a clean arena actions and avoids the need for engineered dissipative
to test thermalization via the Eigenstate Thermalization dynamics—making it uniquely suited to both digital and
Hypothesis [2], and enables direct estimation of free ener- analogue quantum simulation devices. This is achieved
gies and observables that organize phase diagrams [3, 4]. in terms of unitary dynamics, such that an initial ther-
Practically, accurate thermal state preparation on quan- mal state, that can be prepared with reasonable effort,
tum simulators enables us to probe regimes that remain evolves towards a thermal state with respect to a target
challenging for classical methods (e.g., where quantum Hamiltonian. Importantly, this Hamiltonian may well in-
Monte Carlo faces a sign problem [5, 6]). clude interactions that are not native to the system used
Common approaches to preparing thermal states as for implementation. The corresponding unitary dynam-
steady-states of dissipative evolution include using effec- ics will be found as the solution of an optimal control
tive Lindblad dynamics [4, 7–14], or coupling the sys- problem that is independent of the temperature of the
tem weakly to controlled bath degrees of freedom [15– thermal state. As such, the present approach applies to
19], see [20] for an overview. While these protocols have the full range of temperatures.
become increasingly mature, many open questions, such We exemplify our protocol on the paradigmatic cluster-
as the required protocol durations to reach a desired fi- Ising model: an iconic model with a three-body inter-
delity, remain. Yet, there is growing evidence that ther- action that does not exist as a natural interaction in
mal states at low temperature or at criticality require any quantum device. Cluster states, which arise as the
particularly long times to converge to the steady-state, ground states of the cluster–Ising model, appear at fi-
whereas convergence is often much faster at sufficiently nite temperature in contexts ranging from measurement-
high temperature, or far from quantum phase transitions based quantum computing to symmetry-protected topo-
(non-criticality) [21, 22]. logical (SPT) order, and their thermalization properties
These preparation protocols focus predominantly on remain an active area of research [26, 27]. The SPT phase
thermal state preparation for early fault-tolerant quan- and the Ising-ordered phase are separated by a topologi-
tum computers; this requires a complete gate set and cal quantum phase transition [28]. This model is thus ide-
ancillary degrees of freedom, which makes them inappro- ally suited to demonstrate how the present method pro-
priate for analog quantum simulators. Existing protocols vides the theoretical foundation for thermal state prepa-
for these devices focus on thermal state preparation with ration both at criticality and at low temperatures, i.e.,
respect to native interactions: this sharply contrasts with under conditions that proved to be challenging in ear-
analog quantum simulators primarily targetting physics lier approaches. Notably, we observe that the duration
 2

step (i) optimize step (ii) prepare ϱβ (0)

(a) (b)
|z⟩ ∼ p(z)
min
{H(t), K0 } classical sample: for high and low temperature:
1
J
... 0.99

critical

fidelity
|0⟩ |0⟩ |1⟩ |0⟩ = |z⟩

...
pulses or digital: = ϱβ (0) fidelity of ϱβ (tf )
for fixed duration tf

model parameters
step (iii) apply protocol H(t)

optimized
Z1 Z2 Zn−1 Zn

sequence critical

temperature
X1 Xn
from (i) region
...
disordered ordered
X1 X2
step (iv) measure phase phase
Xn−1 Xn

= U |z⟩
...
model parameters
repeat (ii),(iii),(iv)

FIG. 1. The left panel (a) summarises the target thermal state ϱβT preparation as a four-step process. (i) An optimization,
with respect to a figure of merit J , to determine the specific time-dependencies of the system Hamiltonian, H(t), and the
initial condition K0 (Eq. (8)). (ii) The preparation of the thermal state of the initial condition K0 , ϱβ (0), either by sampling
computation basis states as a chain of excited/unexcited spins, |z⟩, from the corresponding Boltzmann distribution, p(z),
(suited to analogue devices), or by a gate sequence (suited to digital devices). (iii) The optimized control sequence, realized
by dynamics induced by H(t) found in (i), is applied to the spin-chain, or implemented via gates in the case of digital devices.
(iv) An observable of interest is measured. Steps (ii),(iii),(iv) are repeated as many times as are needed to accurately capture
the statistics.
The right panel (b) shows a sketch of what can be achieved using the present method. The results in Sec. III C demonstrate that
one can prepare the unitarily evolved thermal state ϱβ (tf ) ≈ ϱβT for a broad range of temperatures, even in critical regions of a
phase diagram. Strikingly, for a given number of spins, the required evolution time tf (to reach a state within some accuracy of
the target) appears to be unaffected by criticality. As such, the results show that one can prepare high-fidelity thermal states
(typically ∼ 0.99 in the worst case) even in critical regimes and low temperature.

of system dynamics required for thermal state prepara- sity matrix

tion is also largely independent of aspects like critical-
ity, and is mostly limited by the fundamental quantum exp(−βKT )
ϱβT = , (1)
speed limit set by the intrinsic system interaction. This tr exp(−βKT )
stands in sharp contrast to dissipative evolution proto-
cols, for which the required evolution time generally in- for the target Hamiltonian KT at inverse temperature
creases substantially in the vicinity of a critical point. β, evolved from an experimentally realizable initial state
Our protocol provides a reliable thermal state prepara- (cf. Sec. III D). Respecting present-day quantum simula-
tion protocol ideally suited to the constraints of present- tion capabilities, we focus on unitary-only protocols gen-
day quantum simulators—relying only on unitary dy- erated by a system Hamiltonian H(t) constrained to ex-
namics generated by native interactions. perimentally available controls—below we consider con-
stant nearest-neighbour interactions and tunable time-
dependent single-qubit terms. In particular, the target
Hamiltonian KT may not lie within the set of available
system Hamiltonians.
II. A BIRD’S EYE VIEW ON THERMAL STATE
In addition, we consider starting from an easy-to-
PREPARATION
initialize thermal state ρβ (0) ∝ exp(−βK0 )—with ini-
tial parent Hamiltonian K0 . Note that the initial state is
Let us first give an overview of our thermal state prepa- prepared at the same inverse temperature β as the target
ration protocol. Our goal is to prepare the thermal den- state, since unitary evolution preserves thermal occupa-
 3

tion, While the restriction to a small Lie algebra enables the

explicit design of system Hamiltonians with optimized
ϱβ (t) = U (t) ϱβ (0) U † (t) ∝ e−βK(t) , (2) time dependence for system sizes that are far out of reach
for descriptions in terms of state vectors or density matri-
i.e., only the parent Hamiltonian K(t)changes under the ces, it poses a fundamental difficulty for optimal control:
any pair of isospectral Hermitian operators are related by
 R
evolution U (t) = T exp −i 0 dsH(s) , with K(0) = K0
t

and a general unitary operator; however, the propagators in-

duced by a Hamiltonian comprised of elements of a small
K(t) = U (t) K0 U † (t) . (3) Lie algebra is necessarily restricted to the corresponding
Lie group. It is thus crucial to choose the initial operator
Thus, to prepare the target state ϱβT = ϱβ (tf ) at final time K(0) = K0 in Eq. (3) such that the target Hamiltonian
tf using unitary evolution, we have to find an experimen- KT can be reached under the restricted system dynamics.
tal protocol H(t) that propagates a suitable initial con- To this end, it is useful to consider a Cartan decom-
dition K0 to the desired target Hamiltonian K(tf ) = KT , position b = k ⊕ m [31, 32] of the Lie algebra b into two
following the von Neumann equation, cf. [29], groups of operators k and m, such that
[k, k] ⊂ k, [m, m] ⊂ k, [k, m] ⊂ m , (6)
K̇(t) = −i[H(t), K(t)]. (4)
where [a, b] ⊂ c is a short-hand notation stating that
The numerical integration of the von Neumann equa- [a, b] ∈ c for any a ∈ a and b ∈ b.
tion (Eq. (4)) is typically prohibited by the exponential Within m, there is a subset h that is maximal Abelian,
growth of a system’s Hilbert space with the number of i.e. all elements of h commute with each other, and there
interacting degrees of freedom. is no further element of m that commutes with all el-
A prominent exception to this rule is given when the ements of h. There is an existence theorem [33] that
system Hamiltonian is comprised exclusively of elements states that for any element m ∈ m there is an element
of a Lie algebra b, i.e. a set of operators {bj } with closed h ∈ h and an element k ∈ k such that
commutation relations
m = exp(ik)h exp(−ik) . (7)
d
[bj , bk ] = −i (5) With an initial condition K0 defined by h and a target
X
λklj bl ,
l=1
KT defined by m (or vice versa), this theorem provides a
guarantee of reachability. As such, there is a guaranteed
in terms of scalar structure constants λklj , where the num- solution to Eq. (3) with K(0) = K0 and K(tf ) = KT
ber of elements (i.e. d) grows sub-exponentially with the at a final time tf , with a unitary induced by a system
number of interacting degrees of freedom, and, is thus Hamiltonian H(t) that generates the Lie algebra b.
smaller than the dimension of the underlying Hilbert Crucially, the initial condition K0 = h is a sum of mu-
space. If the initial condition K(0) can be expanded tually commuting operators, which is of great help for
in terms of elements of the Lie algebra, the solution of the preparation of the initial state ϱβ (0). The existence
Eq. (4), i.e. K(t) in Eq. (3), can be constructed by in- theorem on its own, however, does not specify the initial
tegrating d coupled differential equations, even though condition K0 beyond the restriction to the set h of opera-
construction of the propagator U (t) or the explicit mul- tors. The initial condition from which the control target
tiplication of K(0) and U (t) in Eq. (3) requires an expo- KT is reachable, thus needs to be determined as part of
nentially large effort [30]. the optimal control problem, that can be formalised as
Even if the dynamics U (t)K(0)U † (t) takes place in a min J K(tf ), KT


low-dimensional space spanned by the elements of the H(t), K0 ∈h
Lie algebra, the dynamics U (t)ϱβ (0)U † (t) is typically not s.t. K̇(t) = − i [H(t), K(t)] , (8)
restricted to this low-dimensional space, because the ini-
K(0) = K0
tial thermal state ϱβ (0) includes powers of K0 , which are
not necessarily elements of the Lie algebra. The abil- with an objective function J that is minimal if and only
ity to integrate Eq. (4) efficiently with the initial con- if the objective is reached, and with the time-dependent
dition K(0) = K0 , thus does not typically imply that system Hamiltonian H(t) restricted to practically realiz-
Eq. (4) can be integrated efficiently with an initial con- able interactions and driving protocols.
dition K(0) ∝ exp(−βK0 ). Even though Eq. (4) with The solution of this optimization problem yields the
the initial condition K(0) = K0 can be used to find a desired time-dependent system Hamiltonian H(t) and
time-dependent Hamiltonian H(t) that realises the de- constant parent Hamiltonian K0 . The experimental im-
sired state preparation, the time evolution with a gen- plementation of the dynamics following Eq. (2) that re-
eral initial condition – including the initial thermal states alizes the desired thermal state, then only requires the
ϱβ (0) – can not be simulated classically for large enough ability to prepare the initial state ϱβ (0).
systems, but it needs to be implemented on an actual The complete approach can be summarized in the four
quantum device. steps depicted in Fig. 1a:
 4

(i) an optimization that yields an initial condition K0

and a control sequence that reaches KT ;
(ii) preparation of the initial thermal state ϱβ (0) ∝
exp(−βK0 );
(iii) application of the optimized control sequence;
(iv) measurement of an observable of interest.
Sec. III discusses this general approach more explicitly
with the example of the cluster Ising model.

III. THERMAL STATE PREPARATION OF THE

CLUSTER ISING MODEL

The following discussion will exemplify the general ap- FIG. 2. Schematic representation of the phase diagram of
proach sketched in Sec. II for the specific case of the the cluster Ising model (Eqs. (9), (10)), spanned by the vec-
cluster Ising model based on a system Hamiltonian re- tors ⃗λ = (λ, 0, 0) (left bottom), ⃗λ = (0, λ, 0) (right bottom)
stricted to an extended Ising model. The system Hamil- and ⃗λ = (0, 0, λ) (center top) that define the limiting cases
tonian and control target are defined in more detail in of paramagnetic, Ising-ordered and symmetry-protected topo-
Sec. III A. The Lie algebra generated by the individual logical (SPT) phase.
terms in the system Hamiltonian and its Cartan decom- The model exhibits non-critical behaviour with one dominant
position, as required for reachability and determination term in the Hamiltonian, i.e. domains close to the corners
of the initial condition K0 , is discussed in Sec. III B. Ex- of the triangles depicted in green. If two or all three terms
plicit solutions of optimal control problems are discussed in the Hamiltonian are of comparable magnitude (as indi-
cated in grey), the system shows critical behaviour. In the
in Sec. III C, and protocols for the preparation of the ini-
asymptotic limit n → ∞, the critical and non-critical regions
tial states ϱβ (0) determined by the initial condition are are separated by a sharp phase transition; for finite system
given in Sec. III D. size the crossover between critical and non-critical behaviour
takes place over a finite interval in the parameter space, as
indicated by the shaded domains.
A. System Hamiltonian and cluster Ising model Seven selected points in the parameter space for which opti-
mal control solutions are discussed in Sec. III C are indicated
An explicit implementation of the framework sketched by blue points P1 to P7 , with P1 to P4 in critical and P5 to
in Sec. II requires definition of a system Hamiltonian P7 in non-critical regimes.
H(t) and a target Hamiltonian KT defining the thermal
states to be prepared. of the first and last pair in the chain that resembles the
remains of the three-body XZX interactions obtained
from truncating such a spin chain.
1. Cluster Ising model
The properties of the system eigenstates depend on
the actual values of the scalar constants λ1 , λ2 and λ3
To demonstrate the method’s ability to reliably pre- as sketched in Fig. 2. Because an overall prefactor of the
pare thermal states in critical regimes of Hamiltonians Hamiltonian defines only a natural energy scale, but does
including non-native interactions, the target Hamiltonian not change the physics of the model in any other way,
considered in the following discussion is the cluster Ising the three system parameters can be chosen to adopt the
Hamiltonian normalization
λj = λ , (10)
n n−1
X
KT =λ1 Zj + λ 2 (9)
X X
Xj Xj+1 j
j=1 j=1
n−1 admitting the two-dimensional representation shown in
Fig. 2.
 
− λ3 Z1 X2 + Xj−1 Zj Xj+1 + Xn−1 Zn ,
X

j=2 For λ1 ≫ λ2 , λ3 , the single-qubit Z terms dominate

and the system is deep in a trivial paramagnetic phase.
of an open chain of n spins, comprised of single-spin Z For λ2 ≫ λ1 , λ3 , the two-body XX interactions dom-
terms of magnitude λ1 , nearest neigbour XX interactions inate, and the system is deep in a symmetry-broken
of magnitude λ2 and a three-qubit XZX interaction of Ising-ordered phase. And finally, for λ3 ≫ λ1 , λ2 , the
magnitude λ3 . Because of the open boundary conditions, three-body XZX interactions (together with the two-
there is a Z1 X2 interaction and an Xn−1 Zn interaction body boundary terms) dominate and the system is deep
 5

elements of Lie algebra set index range number of elements row

Zj h⊂m 1≤j≤n n (i)
Xj Zj+1 . . . Zk−1 Xk m 1≤j<k≤n n(n − 1)/2 (ii)
Xj Zj+1 . . . Zk−1 Yk k 1≤j<k≤n n(n − 1)/2 (iii)
Yj Zj+1 . . . Zk−1 Xk k 1≤j<k≤n n(n − 1)/2 (iv)
Yj Zj+1 . . . Zk−1 Yk m 1≤j<k≤n n(n − 1)/2 (v)
Z1 . . . Zj−1 Xj m 1≤j≤n n (vi)
Z1 . . . Zj−1 Yj k 1≤j≤n n (vii)
Xj Zj+1 . . . Zn m 1≤j≤n n (viii)
Yj Zj+1 . . . Zn k 1≤j≤n n (ix)
Z = Z1 . . . Zn h⊂m 1 (x)

TABLE I. Elements of the Lie algebra generated by the operators in Eq. (13) (first column) and sets to which the operators
belong (second column); any element of h is also element of m. The third column depicts the range of the indices of the
operators in the first column, and the fourth column depicts the resultant number of operators.
With n elements each in rows (i), (vi) and (viii), n(n − 1)/2 elements each in rows (ii) and (v) and one element in row (x),
there are n2 + 2n + 1 elements in m. With n(n − 1)/2 elements each in rows (iii) and (iv), and n elements each in rows (vii)
and (ix), there are n2 + n elements in m. The full Lie algebra thus has 2n2 + 3n + 1 elements.
The maximal abelian subset h of m (rows (i) and (x)) has n + 1 elements.

in a symmetry-protected-topological (SPT) phase [28]. spin rotations; it also determines a minimal time required
If two of constants λj are of comparable magnitude, to achieve the thermal state preparation, known as the
the system shows critical behaviour associated with a quantum speed limit [34, 35]. Modulating the single
quantum phase transition. The energy gap between the qubit energies Zj enables effectively weakening of the
ground state and first excited state is finite for any fi- XX interactions, similar to the process of dynamical de-
nite system size, but the gap closes asymptotically with coupling. The three-body XZX interactions in Eq. (9)
a power-law in the large-system limit. The critical and can be obtained with this system Hamiltonian as effective
non-critical regions are depicted in the triangular phase third-order processes owing to the commutation relation
diagram (Fig. 2), where each corner indicates a dominant
phase. Because the locations of quantum phase transi- [Xj−1 Xj , [Xj Xj+1 , Zj ]] = 4Xj−1 Zj Xj+1 . (12)
tions are not sharply defined, we shade a broad region
to indicate their approximate positions. The additional The system Hamiltonian H(t) (Eq. (11)) thus seems
features of Fig. 2 are discussed in Sec. III C. sufficiently elementary to be implemented with suitable
time-dependence in practice, but also sufficiently general
to induce the dynamics required to emulate the cluster
2. System Hamiltonian Ising model of Eq. (9).

While the single qubit terms Zj and the interactions

Xj Xj+1 can be realized directly in many platforms of B. Lie algebra and Cartan decomposition
well-controllable qubit registers, the three-body interac-
tions Xj−1 Zj Xj+1 generically needs to be realized via The Lie algebra generated by the operators
effective processes, and thus it cannot be assumed intrin-
sic to realistic hardware and should not be included in X1 , (13a)
the system Hamiltonian. Xn , (13b)
In the following, the system Hamiltonian is thus con- Zj , (13c)
sidered to be of the form
Xj Xj+1 , (13d)
n−1 n
(x)
H(t) = g Xj Xj+1 + hj (t)Zj + hj (t)Xj , contained in the system Hamiltonian H(t) (Eq. (11))
X X X

j=1 j=1 j∈{1,n} is a maximal set of linearly independent operators that

(11) includes all terms in Eq. (13), commutators and nested
with a nearest neighbour XX interaction of constant commutators of these terms. For example, the operator
magnitude g, tuneable time-dependent single-spin ener- Y1 X2 is an element of the Lie algebra, since it is obtained
gies Zj and additional tuneable X driving of the end via the commutator of the single qubit term Z1 and the
spins in the chain. The interaction constant g defines a two-qubit interaction X1 X2 . The interaction term Z1 Z2 ,
natural timescale of any system dynamics beyond single- on the other hand is not part of the Lie algebra, since it
 6

can not be obtained in terms of nested commutators of

the terms in Eq. (13).
A Cartan decomposition can be defined with the help
of an involution Λ, i.e., a linear map such that Λ2 is
the identity map and such that Λ([A, B]) = [Λ(A), Λ(B)]
for any pair of operators A and B. The commutation
relations Eq. (6) are necessarily satisfied, if the elements
k ∈ k satisfy the eigenvalue relation Λ(k) = k and the
elements m ∈ m satisfy the eigenvalue relation Λ(m) =
−m.
For the specific Lie algebra generated by Eq. (13), each
element is a string of Pauli operators, and a Cartan de-
composition can be defined in terms of the transposition
FIG. 3. Thermal state infidelity (Eq. (19)) as a function of
Λ(A) = −AT . Since the operator Y satisfies the relation
λβ, with the scale λ given in Eq. (10), for the two exemplary
Y = −Y T , while X and Z satisfy the relations X T = X points P1 and P7 depicted in the phase diagram in Fig. 2
and Z T = Z, the operators k have an odd number of fac- for a system of n = 11 spins. The optimal control solution
tors Y , whereas the operators m have an even number of for the thermal state preparation is obtained by optimizing
factors Y . for the operator infidelity, and the final value obtained in the
Table I depicts all the elements of the Lie algebra gen- numerical optimization is depicted as a dashed line.
erated by Eq. (13) and specifies to which set in the Cartan The point P1 (left inset) corresponds to critical behaviour
decomposion each element belongs. A maximal Abelian with a closing gap; the state fidelity at low temperatures is
subset h ⊂ m is given by thus larger than the operator fidelity in contrast to the case
of non-critical behaviour depicted in the right inset. In both
h = {Z1 , . . . , Zn , Z} (14) cases, it is evident that the actual figure of merit (the thermal
state infidelity J (ϱβ (tf ), ϱβT )) remains small even though the
including all the single-qubit operators Zj and the n- proxy (the operator infidelity J (K(tf ), KT )) was optimized.
body interaction
n
Z=
Y
Zj , (15) With the cluster Ising model defined in Eq. (9), the
system Hamiltonian defined in Eq. (11), and the initial
j=1
condition K(0) = K0 restricted to the set h in Eq. (14),
that will be referred to as parity in the following. the general optimization problem given in Eq. (8) is com-
For the optimization (Eq. (8)) this means that the ini- pletely specified with the definition of the objective func-
tial condition K0 is of the form tion J and the parametrization of the time-dependencies
n
of the functions hj (t) to be optimized.
K0 = c0 Z +
X
cj Zj (16) Because of its practical evaluability, all of the explicit
results discussed in the following are based on the infi-
j=1
delity
with scalars cj to be optimized over.
tr AB
The initial state ϱβ (0) for an experimental implementa- J (A, B) = 1 − p (17)
tion of the thermal state ∝ exp(−βKT ) following Eq. (2) tr (A2 ) tr (B 2 )
is given by ϱβ (0) ∝ exp(−βK0 ) with the values of the as the explicit choice for the objective function
coefficients cj obtained from the solution of the optimiza- J (K(tf ), KT ) in Eq. (8).
tion (Eq. (8)). A brief discussion of subtleties resultant from this
Moreover, the initial states ϱβ (0) to be prepared are choice of objective function, the explicit parametrization
generally diagonal in the computational basis defined of time-dependencies as required for the numerical imple-
by the operators Zj , but the n-body interaction term mentation of the optimization, and the implementation
Z implies that the initial thermal states are generally used for this work are provided in App. B.
not product states. Nevertheless, the states can be pre- Explicit control targets are thermal states correspond-
pared rather efficiently as discussed in more detail in ing to parameters depicted in Fig. 2 with points P5 , P6 ,
Sec. (III D). and P7 in a regime with one single dominant term in the
target Hamiltonian; points P2 , P3 , and P4 in a regime in
which two terms are of comparable magnitude and the
C. Control results
third term of much smaller magnitude; and the point P1
in a regime with all the three terms of the Hamiltonian
In the following, we outline the optimization results of comparable magnitude.
depicted by step (i) in Fig. 1, and for moderate system A vanishing value of the operator infidelity
sizes, we explicitly simulate the thermal state dynamics
depicted by step (iii) in Fig. 1. J (K(tf ), KT ) (18)
 7

(with λ defined in Eq. (10)) for the exemplary cases P1

and P7 . A dashed line depicts the value of the operator
infidelity obtained for this optimization. The exemplary
cases P1 and P7 include one critical and one non-critical
point in the phase diagram Fig. 2. Corresponding data
for P2 to P6 is shown in Fig. 8 in App. A to provide
evidence for qualitatively similar behaviour of all critical
points and of all non-critical points.

The ground state infidelity (i.e. the limit β → ∞) is

larger than the operator infidelity for P1 , but it is smaller
than the operator infidelity for P7 . This can be attributed
to the fact that the control target KT has a sizeable
gap for P7 , whereas the gap for P1 is closing as a power
law in the total number of spins. Since the accuracy
of the ground state fidelity is bounded by the operator
fidelity with a bound that depends on the gap [30] it is
expected that control targets with a small gap require
lower operator infidelities to reach a given ground state
infidelity than targets with a larger gap.

At sufficiently high temperatures, the thermal state

FIG. 4. Operator infidelity resultant from the optimization infidelities decrease as the thermal state approaches full
for the two exemplary points P1 and P7 depicted in Fig. 2 degeneracy; here, imperfections in the state transfer of
as a function of system size (green diamonds). State fideli- individual eigenstates are increasingly averaged out by
ties obtained with the corresponding optimal control solutions the nearly uniform populations. At in-between tempera-
are depicted for system sizes up to n = 11 spins both for the tures, a sizeable number of excited states contributes to
ground state (β → ∞, red triangles) and for the finite temper-
the thermal state infidelity, but the thermal states are
ature β = 2/λ (purple crosses) with the scale λ (Eq. (10)) of
KT . The ground state infidelity upper bound Bgs (Eq. (20))
not close to degeneracies. The infidelities for such ther-
is depicted by empty square points. mal states can exceed ground state infidelity, and this
Similar to Fig. 3, state infidelities tend to be larger than op- is typically the case for states in non-critical regions in
erator infidelities in critical regions (top panel) and smaller which low ground state infidelities are obtained.
in non-critical regions (bottom panel). In both panels, the
operator infidelity is a sufficiently good proxy for the state Since the comparison between operator infidelity and
infidelity to be used as control target. thermal state infidelity depends also on the system size,
Fig. 4 depicts infidelities achieved with optimized proto-
cols for the points P1 and P7 as a function of system size
verifies that the final operator K(tf ) matches the control n (corresponding data for P2 to P6 is shown in Fig. 9 in
target KT exactly. Even though an exactly vanishing in- App. A). Thermal state infidelities with the intermediate
fidelity indicates that the actual goal – realization of the temperature β = 2/λ are depicted for system sizes up
desired thermal state – can be achieved perfectly, mini- to n = 11. As seen in Fig. 3, most prominently for the
mization of the operator infidelity is not the actual goal non-critical point, this temperature typically represents
of the optimal control problem; it is rather a substitute a worst-case regime for the thermal state infidelity. Even
problem that is advantageous to implement, but the ac- in this regime, the thermal state infidelities are typically
tual figure of merit of interest is the state infidelity within an order of magnitude higher than the correspond-
ing operator infidelity for targets in critical regimes. For
J (ϱβ (tf ), ϱβT ) (19) non-critical targets, the thermal state infidelity is gen-
erally within an order of magnitude lower than the cor-
of the system state ϱβ (tf ) (Eq. (2)) at the final time tf responding operator infidelity. For targets in the critical
with respect to the target state ϱβT (Eq. (1)). Given the region, the ground state infidelity – i.e. the thermal state
exponential scaling of the Hilbert space, this infidelity infidelity in the limit β → ∞ – commonly represents
can be evaluated in practice only for systems with a mod- the worst-case thermal state infidelity due to the clos-
erate number of spins, such as n = 11, as is the case in ing ground state energy gap. This is particularly evident
Fig. 3. However, owing to the polynomial scaling of the in Fig. 9 in App. A. Like the intermediate-temperature
Lie algebra (Tab. I), Eq. (18) can be optimized for sys- state fidelity, this quantity is infeasible to evaluate for
tem sizes beyond what is classically simulable in terms large system sizes (n ≥ 14). Therefore, it is instructive
of thermal state dynamics. In such cases quantum hard- to introduce an upper bound for the ground state infi-
ware is required to explicitly verify the state. delity that can be evaluated even when explicit classical
Fig. 3 depicts the state fidelity as a function of λβ simulation is infeasible. Proven in [30], the upper bound
 8

is given by ulate the system dynamics, the effort to optimize for

the time-dependent system Hamiltonian H(t) grows with
⟨Ψ(tf )| KT |Ψ(tf )⟩ − E0 the system size n (i.e., Eq. (8) requires the optimiza-
Bgs = , (20)
E1 − E0 tion of O(n2 ) parameters). In contrast to the simula-
tion of the system dynamics that mostly requires com-
where E0 and E1 are the ground state and first ex- putational memory (spatial complexity), the optimiza-
cited state energies of KT respectively, and |Ψ(tf )⟩ is the
tion requires the ability to try many different time-
ground state of the evolved operator K(tf ). The expecta-
dependencies; since this can be done sequentially, this
tion value ⟨Ψ(tf )| KT |Ψ(tf )⟩ can be evaluated implicitly,
implies time-complexity, i.e., a cost dominated by the
without explicit construction of KT and the full statevec-
total runtime needed to explore many candidate time-
tor |Ψ(tf )⟩ [30].
dependencies rather than by memory usage. As such,
Strikingly, in Fig. 4, it appears that the bound tightens
one can always attempt to find a better time-dependence
with respect to the operator infidelity as the system size
than the best found so far, but in practice one will need
increases. For the critical example, P1 , it appears that
to accept a solution with some finite accuracy. As seen
the bound is substantially lower for n = 17, and espe-
in Figs. 3 and 4, operator infidelities of O(10−4 ) can be
cially n = 14. However, this indicates the existence of an sufficient to design practical protocols for thermal state
optimal control solution with particularly good ground preparation, but an unambiguous identification of the
state infidelity. quantum speed limit requires substantially lower infideli-
Crucially, Figs. 3 and 4 demonstrate that thermal ties. The numerical assessment of this time-scale is there-
states in the critical region can be prepared with high- fore limited to at most n = 10 spins.
fidelity using the same experimental resources–like evo-
For the practically accessible system sizes, however,
lution time tf –required for non-critical regions.
Fig. 5 (and Fig. 10 in App. A) allows us to read off
All the infidelities shown in Figs. 3 and 4 are finite,
the system-size dependent quantum speed limit rather
but sufficiently small to realize thermal state prepara-
accurately. This duration is depicted for all the points
tion with accuracies compatible with existing or near-
P1 to P7 in Fig. 6 as a function of the system size n.
term hardware. Fundamentally, the achievable value of
Strikingly, the differences of the minimal duration ob-
infidelity is limited by the duration tf of the controlled dy-
tained for different points Pj at a given system size are
namics, because the interaction constant g in the system
within the accuracy to which this time can be read off
Hamiltonian H(t) (Eq. 11) defines a fundamental mini-
of Figs. 5 and 10. There is thus no significantly longer
mal time scale (quantum speed limit) for any entangling
duration required for the preparation of states in critical
dynamics.
regimes than for non-critical ones. This contradicts the
While optimal control solutions can overcome restric-
common perception that the preparation of critical states
tions of other existing strategies, such as adiabaticity, the
is limited by the time-scale imposed by the closing gap.
restriction imposed by the quantum speed limit is an in-
This perception, however, is based on adiabatic dynam-
trinsic system property that can not be overcome simply
ics. The results in Fig. 6 are thus not in contradiction
in terms of sufficiently strong or sophisticated driving of
with the common perception, but they are an indication
single-qubit dynamics.
of the non-adiabatic character of the control solutions
As such, one would expect that there is a minimal
that can be found with the present approach.
duration tmin for the evolution time tf that is required
To the extent that Fig. 6 allows for an identification
for preparation protocols with vanishing infidelities or,
of a scaling relation, the minimal time required for the
in practice, values of infidelity that are limited only by
state preparation appears to grow linearly with the sys-
the accuracy of the numerical opimization.
tem size. This is consistent with the scaling obtained for
Fig. 5 depicts values of the operator infidelity as a func-
the preparation of the ground state of the cluster Ising
tion of the evolution time for various system sizes for the
Hamiltonian with λ1 = λ2 = 0 in Eq. (9) [30]. This
critical points P1 to P4 (Fig. 2). (Data for the noncritical
scaling also highlights that the state preparation is not
points P5 to P7 is shown in Fig. 10 in App. A). Apart
relying on an adiabatic approximation, that would yield
from numerical noise caused by the common limited reli-
a scaling with the closing gap.
ability of numerical optimizations, the operator infidelity
decays monotonically with increasing evolution time tf .
This decay is rather weak for sufficiently short times, but
there is a narrow time-window, in which we observe a D. Explicit initial thermal state construction
very pronounced drop of infidelity, down to a value that
is limited by the numerical accuracy of the optimization. In addition to the explicit control solutions constructed
The clearly pronounced sudden drop in operator in- to assess the infidelities discussed in Sec. III C above, an
fidelity suggests that the optimization is indeed able to implementation on a quantum device also needs a proce-
find close-to-perfect control solutions, as soon as the evo- dure to prepare the initial thermal state ϱβ (0), as seen in
lution time tf is sufficiently long to admit such solutions. Eq. (2). This part of the broader method is depicted as
While the present approach allows us to overcome step (ii) in Fig. 1.
the exponential growth of the numerical effort to sim- The construction of optimal control solutions in
 9

FIG. 6. The quantum speed limit duration, gtmin , where g

FIG. 5. Best operator infidelity obtained as result of the op- is the interaction constant in the system Hamiltonian H(t)
timization as a function of the evolution time gtf , where g (Eq. (11)), as a function of system size n. tmin is identified
is the interaction constant in the system Hamiltonian H(t) as the clear drop in operator infidelity shown in Fig. 5 and
(Eq. (11)). There is a well-defined drop in the infidelity, sug- Fig. 10 in the supplementary solutions. The data suggests a
gesting that there is a minimal duration of system dynamics linear relationship between tmin and system size, regardless of
required to realize the thermal state preparation. This mini- phase diagram location (Fig. 2), which is a clear indication of
mal duration appears independent of the specific target within non-adiabaticity. The phase diagram (Fig. 2) is included as a
the cluster Ising model, and it seems to be growing only mod- visual aid, along with a dashed linear reference line.
erately with the system size n.

1. The initial state

As discussed in Sec. III B, the initial condition K0 is

of the form
Sec. III C is in the spirit of an analog quantum simu- n
lation with a Hamiltonian (Eq. (11)) that has a constant K0 = c0 Z + (21)
X
cj Zj
interaction and additional single-spin driving that is be- j=1
ing applied while the spins are interacting. The frame-
work of optimization is, however, also directly applica- with the weights ci of the individual terms determined
ble to a digital quantum simulation with two-qubit XX- as solutions of the optimal control problem. The initial
gates and single-qubit gates that are being applied se- condition for an actual implementation is then a thermal
quentially. While the present framework to optimize the state
thermal state preparation is applicable to both analog
and digital quantum simulations, the initial state prepa- exp(−βK0 )
ϱβ (0) = . (22)
ration is likely to be different in the two approaches. The tr exp(−βK0 )
following discussion of initial thermal state preparation
thus includes two different approaches. with the initial condition K0 (Eq. (21)) as the parent
Hamiltonian.
Without the parity Z, i.e. the n-body interaction, the
The first approach realizes the desired initial thermal initial condition ϱβ (0) would be a product of single-qubit
states in terms of an average over randomly selected pure thermal states that is rather straightforward to prepare.
states. This approach does not require the application The following discussion will thus mostly focus on the
of any entangling gates, and is thus suitable for analog aspects required to prepare a state that is not a product
quantum simulations. The second approach is defined state.
in terms of a gate sequence that realizes a pure state
with the occupations that are consistent with the desired
mixed intial state, and that turns this coherent superpo- 2. Random sampling
sition into an incoherent mixture with additional entan-
gling gates that involve additional qubits that play the Any mixed state can be realized in terms of an ensem-
role of an environment. This approach does require the ble of pure states with associated probabilities. Such
capabilities of a digital quantum simulator. probabilities can easily be constructed for the initial
 10

ϱβ (0) (Eq. (22)), but in order to realize the state prepa- |Ψ⟩
ration in an experiment, it is important that one can
sample efficiently from the exponentially large ensemble. Ry (2ϑ1 )
A precondition for finding such a sampling procedure is .. .. ..
S ϱβ (0)
the construction of the probabilities associated with each . . .
pure state in the ensemble. It is natural to take these Ry (2ϑn )
pure states to be the eigenstates of the initial operator
K0 (Eq. (21)), i.e. the product states
A .. .. trA
|⃗zn ⟩ = |z1 ⟩ ⊗ . . . ⊗ |zn ⟩ (23) . .

of single-spin Zj eigenstates. Uχ
M
For any traceless Hamiltonian with eigenvalues ±c, the
probability p(z) for a given pure state with eigenvalue zc
and z ∈ {−1, 1} in an ensemble representing a thermal FIG. 7. Quantum circuit diagram to realize the state in
state is given by Eq. (33). S and A are n qubit registers, representing sys-
tem and auxiliary Hilbert spaces, and M is a single qubit to
1 exp(−βcz) encode parity. All registers are initialised in |0⟩. The gate
p(z) = , (24)
nd exp(−βc) + exp(βc) sequence in the yellow box generates the pure state |Ψ⟩, after
which discarding the auxiliary system yields, trA |Ψ⟩ ⟨Ψ| = ϱβF
where nd is the degree of degeneracy of the eigenvalues. (Eq. (35)). A unitary Uχ = |0⟩ ⟨χ| + |1⟩ ⟨χ⊥ | applied on qubit
Given that the quantum number z adopts only the values M , followed by a Z-measurement, yields the desired thermal
±1, this can be written equivalently as state (Eq. (33)) in register S , where {|χ⟩ , |χ⊥ ⟩} is an or-
thonormal basis. The quantum circuit requires n single qubit
1 1 − mz Y rotations (with specific angles given in Eq. (39)), and 2n
p(z) = Q(m, z) , with Q(m, z) := , (25)
nd 2 CNOTs, where the first n CNOTs can be applied in a single
unit of depth.
where m = tanh(βc).
Since the operator K0 (Eq. (21)) is a sum of mutually
commuting terms, the probability for a given state |⃗zn ⟩ of the normalization constant Nβ .
in the ensemble representing the thermal state ϱβ (0) can With the marginal probabilities pj (⃗zj ) (Eq. (27)) for
be written as the states of the first j spins, one also obtains the condi-
tional probabilities
1 Y
n
p(⃗zn ) = Q(mj , zj ) , (26)
Nβ j=0 p(⃗zj+1 )
p(zj+1 |⃗zj ) = (30)
p(⃗zj )
with mi = tanh(βci ) and z0 = j=1 zj ; this is essentially
Qn
the product of the probabilities resultant from each term for the spin state |zj+1 ⟩ given a state |⃗zj ⟩ of the first j
in the Hamiltonian, but since the quantum number z0 spins, and these conditional probabilities give access to
associated with the operator Z is not an independent an efficient sampling procedure.
variable, the normalization Nβ is not just the product of In order to pick a given state |⃗zn ⟩ with the appropri-
the normalizations of the individual probabilities. ate probability, one can pick a state |z1 ⟩ of the first spin
Leaving aside that the normalization factor Nβ is still with the probability p(⃗z1 ) given in Eq. (27). The prob-
undetermined, oneP can construct the marginal probabil- ability for a state |z2 ⟩ of the second spin is given by the
ities p(⃗zj−1 ) = zj =±1 p(⃗zj ), and one obtains conditional probability p(z2 |⃗z1 ) (Eq. (30)), which can be
explicitly evaluated because a state |z1 ⟩ for the first spin
1
j j
is selected. This autoregressive sampling follows through
!
(j)
p(⃗zj ) = Q(mi , zi ) (27)
Y Y
Nβ
Q m0 , zi to the last spin.
i=1 i=1
This procedure allows one to sample pure states |⃗zn ⟩
with from the probabilities in Eq. (26). An average over sam-
ples created in this way realizes the desired thermal state
n
(j) ϱβ (0).
m0 = (−1) (28)
n−j
Y
m0 mi .
i=j+1

The condition p(⃗z1 ) = 1 finally determines the

P
z1 =±1
3. Gate circuit
value
Alternatively, instead of the classical sampling dis-
1
n
!
Nβ = 1 − (−1)
n
(29) cussed above, the initial thermal state ϱβ (0) can be
Y
mi
2 i=0 prepared by realizing system-environment entanglement
 11

with a sequence of gates, akin to the preparation of ther- With tr ϱβNI Z = (−1)n i=1 tanh(βci ) and ⟨χ| Z |χ⟩ =
Qn
mofield double states [36, 37]. Although the initial con- − tanh(βc0 ), this reduces to
dition K0 (Eq. (21)) includes the Z-string Z, so ϱβ (0)
cannot be prepared as the product of single-spin mixed 1
n
!
Ps = 1 − (−1)n tanh(βci ) . (38)
Y
states, it is nevertheless instructive to first follow the 2
preparation of single-spin mixed states. i=0

This can be done with the system-spin prepared in the In the high temperature limit, the success probability
state |0⟩, a Y -rotation by an angle determined by the approaches the value of 1/2. In the low temperature
temperature, and a CNOT gate with the system spin limit, the successPprobability is close to unity if the non-
as control qubit and an auxilary spin as target qubit to interacting part j=1 cj Zj has the correct parity to also
n
realize the system-environment entanglement. With n be a ground state of c0 Z. In the opposite case the success
system spins and n corresponding auxilary spins, this probability can also be too low to be practical. In this,
allows one to realize any thermal P state ϱβNI with a non- one can use the deterministic sampling methods instead.
interacting initial condition K0 = j=1 cj Zj , i.e. c0 = 0
n
Fig. 7 depicts a circuit diagram to realize this state
in Eq. (21). preparation. It consists of n qubits labelled S to realize
In the general case of a finite parity component in the the system spins, n qubits labelled A to serve as auxiliary
initial condition (i.e. c0 ̸= 0), the initial thermal state to spins and one additional qubit labelled M required to
be prepared is of the form realize the Boltzmann factor with the total parity Z.
All system qubits are initialized in their state |0⟩, and
ϱβ (0) ∝ ϱβNI e−βc0 Z . (31) single qubit rotations exp(−iϑj Yj ) with the angles ϑj ,
such that
With the identity
1 − mj 1 + mj
r r
cos ϑj = sin ϑj = , (39)
e−βc0 Z = e−βc0 P+ + eβc0 P− , (32) 2 2

in term of the projectors P± = 12 (1 ± Z) onto positive with mi = tanh(βci ) on each of the qubits S ensure
and negative parity, and commutativity [Z, ϱβNI ] = 0, this that these qubits have the occupations that are consistent
reduces to with the single-spin thermal state ∝ exp(−βcj Zj ).
A subsequent series of CNOT gates with qubits S as
ϱβ (0) ∝ e−βc0 P+ ϱβNI P+ + eβc0 P− ϱβNI P− . (33) control qubits and qubits A as target qubits ensures that
the reduced density matrix of the system of qubits S is
The desired state is thus a weighted, incoherent sum of given by an incoherent mixture and not by a coherent
superposition.
the two parity components P+ ϱβNI P+ and P− ϱβNI P− of the
A further series of CNOT gates with the qubits S as
non-interacting state ϱβNI . control qubits and the qubit M as target qubit is a first
In order to realise such a mixture in a gate circuit, it is step in order to adjust the occupations of the individ-
helpful to introduce an additional spin M , and to apply n ual components in this mixture to take into account the
CNOT gates with the individual system spins as control parity term ∝ exp(−βc0 Z). After a subsequent Y rota-
qubits and the spin M as target. tion exp(iϕY ) with tan(ϕ) = eβc0 on qubit M (equivalent
Starting with ϱβNI ⊗ |0⟩ ⟨0|M , this yields the state to the application of the unitary Uχ depicted in Fig. 7),
projection onto the state |0⟩ upon a readout of qubit M ,
ϱβF = ϱβNI P+ ⊗ |0⟩ ⟨0|M + ϱβNI P− ⊗ |1⟩ ⟨1|M (34) leaves the reduced system of qubits S in the desired ther-
1 β  mal state including the parity contribution.
= ϱNI ⊗ 1M + ϱβNI Z ⊗ ZM . (35)
2
A subsequent measurement on spin M in a basis includ- IV. DISCUSSION AND OUTLOOK
ing the state
βc0 βc0
The present approach enables quantum simulations
e− 2 |0⟩ + e 2 |1⟩ that require both the realization of effective processes
|χ⟩ = (36) beyond native interactions and state preparation beyond
2 cosh(βc0 )
p
the zero-temperature limit. As we demonstrate using the
yields the desired thermal state (Eq. (33)) conditioned cluster Ising model thermal state as the control target,
on the projection of the spin M onto the state |χ⟩. our method works reliably both throughout the entire
The success probability Ps of this protocol is given by phase diagram and across all temperature regimes. Un-
the norm of the projected state, i.e. like previous approaches, our results show that the ther-
mal state evolution duration tf is strikingly unaffected by
1  aspects like criticality or temperature. While the cluster
Ps = tr ⟨χ| ϱF |χ⟩ = 1 + (tr ϱβNI Z) ⟨χ| Z |χ⟩ . (37) Ising model is an example to illustrate the workings of
2
 12

the present framework, the methodology is directly ap-

plicable to any thermal state of a Hamiltonian within the
Lie algebra discussed in this paper or any other Lie alge-
bra with suitable scaling – some examples can be found
in [38].
Aspects like disorder in the system Hamiltonian H(t)
and in the target Hamiltonian KT are within the scope of
a given Lie algebra [39] and can thus be easily taken into
account in the construction of explicit solutions. Also,
the restriction to thermal states in this work is due to
their practical importance, but the framework is directly
applicable also to mixed states that are not of Gibbs-
form.
Any problem that does not fit into the framework of a
small Lie algebra does require numerical approximations
in order to benefit from the methodology laid out here.
Additional terms that exceed the utilized Lie algebra can
be taken into account perturbatively [39] without jeop-
ardizing the numerical efficiency that lies at the core of
the present framework.
While the explicit optimizations presented here are in
the spirit of an analogue quantum simulation, the for-
malism can equally well be applied to the optimization
of gate parameters for a digital quantum simulation. In FIG. 8. Thermal state infidelity (Eq. (19)) as a function of λβ,
such a case, there is no intrinsic interaction constant to with the scale λ given in Eq. (10), for the points P2 , . . . , P5
define a fundamental time scale, and the figure of merit depicted in the phase diagram in Fig. 2 for a system of n = 11
replacing the overall system evolution time could be the spins. The optimal control solution for the thermal state
preparation is obtained by optimizing for the operator infi-
number of entangling gates that need to be applied to a
delity, and the final value obtained in the numerical opti-
pair of qubits before infidelities on the level of numerical mization is depicted as a dashed line.
accuracy can be reached. The observations are largely the same as the exemplary points
Practical considerations such as addressability of indi- (Fig. 3): the actual figure of merit (the thermal state infidelity
vidual qubits or spectral properties of the driving func- J (ϱβ (tf ), ϱβT )) remains small even though the proxy (the op-
tions can be taken into account in the optimization. The erator infidelity J (K(tf ), KT )) was optimized. For critical
explicit control solutions presented here, assume that points P3 , P4 , the thermal state infidelity in the low tempera-
each of the qubits can be driven individually, but if an ture regime can exceed the operator infidelity, however given
actual device cannot achieve perfect addressing, or even a small operator infidelity solution, the worst-case thermal
only global control, then the framework can be applied state infidelity remains sufficiently accurate.
to a system Hamiltonian that respects such a restriction.
With the flexibility to accommodate practical require-
ments and to realize challenging quantum simulations in the Imperial HPC cluster.
a fashion that requires minimal resources from a quan-
tum device, the present approach allows us to push the
envelope of quantum simulation further towards practical DATA AVAILABILITY
use-cases.
The optimal control solutions used in this paper are
available without restriction [40].
ACKNOWLEDGEMENTS

This work was supported by the U.K. Engineering Appendix A: Supplementary control solutions
and Physical Sciences Research Council (EPSRC DTP
- EP/W524323/1). MB and PMS were funded by the In this section, we show the control solutions for the
European Union (ERC, QuSimCtrl, 101113633). Views remaining points of the phase diagram, P2 . . . , P6 . As
and opinions expressed are however those of the authors labelled in Fig. 2, points P2 , P3 , P4 are critical, and points
only and do not necessarily reflect those of the Euro- P5 , P6 are non-critical.
pean Union or the European Research Council Executive While the non-critical points typically have a sizeable
Agency. Neither the European Union nor the granting ground state energy gap that does not close with system
authority can be held responsible for them. Numerical size, the target KT corresponding to parameters speci-
simulations and optimization routines were performed on fied by the non-critical point P6 has a notable spectrum.
 13

FIG. 10. Best operator infidelity obtained as result of the

optimization as a function of the evolution time gtf , where g
is the interaction constant in the system Hamiltonian H(t)
(Eq. (11)), for the non-critical points P5 , P6 , P7 . The location
of the well-defined drop in the infidelity, which numerically
identifies the quantum speed limit, is in close agreement with
the data given Fig. 5.

For odd system sizes n, the gap is thermally realizable;

however, it remains small and constant for all odd n –
approximately an order of magnitude below the gap for
critical targets (i.e., points P1 , . . . , P4 in Fig. 2) at large
system sizes, such as n = 26. As a result, the numeri-
cal deviations of K(tf ) from the true target KT must be
sufficiently small, so that the perturbation does not in-
duce substantial mixing among the thermally populated
FIG. 9. Operator infidelity resultant from the optimization low-lying eigenstates of KT .
for the points P2 , . . . , P6 depicted in Fig. 2 as a function of
system size (green diamonds). State fidelities obtained with To that end, we note that the small gaps for the point
the corresponding optimal control solutions are depicted for P6 lack physical relevance – they arise from a small
system sizes up to n = 11 spins both for the ground state perturbation from the exactly two-fold degenerate pure
(β → ∞, red triangles) and for the finite temperature β = Xj Xj+1 Hamiltonian rather than any form of critical-
2/λ (purple crosses) with the scale λ (Eq. (10)) of KT . The ity. Because the only requirement is finding a more accu-
ground state infidelity upper bound Bgs (Eq. (20)) is depicted rate optimal control solution, rather than increasing the
by empty square points. Limited data is shown for non-critical physical evolution tf , this is a matter of optimization-
point P6 , justified in App. A. Similar observations are made time complexity, and thus the present method can still
to Fig. 4: in all cases, the operator infidelity is a sufficiently
be used to prepare a target thermal state in such regions,
good proxy for the state infidelity to be used as control target.
provided the gap is large enough that the ground state
can still be distinguished numerically.
Nevertheless, because of the unusually small gaps, we
In particular, for even n system sizes, the gap closes to only consider the intermediate temperature state infi-
numerical accuracy even at moderate system sizes. For delity for this point. This point is thus omitted from
example, with n = 14 and the scale λ = 1 (Eq. (10)), Fig. 8 (the supplementary plot of Fig. 3), and limited
the gap is O(10−13 ). At intermediate temperature, the data is shown in Fig. 9 (the supplementary plot of Fig. 4).
ground state is a maximal mixture of the ground and first As with points P1 and P7 in the main text, the operator
excited state to a very good approximation. As such, the infidelity J (K(tf ), KT ) is evidently a good proxy, ensur-
thermal state fidelity can still be very accurate. However, ing that the thermal state infidelity J (ϱβ (tf ), ϱβT ) is also
because of the size of the gap, it is not computationally small, and the plots show qualitatively similar behaviour
feasible to thermally realize the actual ground state and for points in critical and non-critical regions, with P6 as
compute the infidelity for the limit β → ∞. an extreme exception.
 14

Fig. 10 is the supplementary plot of Fig. 5 and shows where a(t) is a time-dependent coefficient vector with
the drops in operator infidelity required to resolve the elements {aj (t)} with length equal to the number of el-
quantum speed limit duration for the remaining non- ements in the Lie algebra, and G(t) is a matrix with
critical points, such that the durations can be read-off elements,
the plot. Strikingly, the specific speed limit durations
tmin appear to coincide with that of critical points, and
Glj (t) = hk (t)λklj . (B9)
X
additionally grow linearly with system size (as shown in
Fig. 6).
k

To that end, the differential equation (Eq. (B8)) provides

Appendix B: Details of numerical implementation an efficient numerical representation of the von Neumann
equation, owing to the quadratic scaling of elements in
Here, we provide additional details for the technicali- b. Additionally, because G(t) is sparse, Krylov methods
ties of the optimization including the equation of motion can be employed for further computational efficiency.
required for efficient numerical propagation using the von
Neumann equation (Eq. (4)) discussed in App. B 1, and
analytical gradients derived in App. B 2. Such details are
largely based on the standard implicit control framework 2. Analytical gradient of the objective function
[30]. Finally, we end this section with some more subtle
technical details including further details of the objective As is common practice with GRAPE-inspired opti-
function and discretization of the system Hamiltonian’s mization, it is often the case that analytical gradients
time-dependence in App. B 3. enable more efficient convergence than numerical gradi-
ent estimators such as finite differences. By defining a
target coefficient vector aT with elements {aT,j } such
1. Equation of motion that the control target can be written using the expan-
sion,
The control terms,
KT = (B10)
X
{Zi , X1 , Xn , Xj Xj+1 } , (B1) aT,j bj ,
j
where i ∈ {1, . . . , n} and j ∈ {1, . . . , n − 1}, form the
polynomially scaling Lie algebra b, with elements {bj } the operator infidelity J (K(tf ), KT ) (Eq. (17)) takes the
satisfying closed commutation relations (Eq. (5)). form,
Writing the system Hamiltonian as the expansion,

H(t) = hj (t)bj , (B2) a(tf ) · aT

X
J =1− , (B11)
j ∥a(0)∥∥aT ∥
and,
where we have used that ∥a(tf )∥ = ∥a(0)∥. Eq. (B11)
K(t) = aj (t)bj , (B3) which can be efficiently evaluated.
X

j We provide two derivations: the analytical gradient

of J with respect to the discretized time-dependence of
the von Neumann equation can be re-written as the system Hamiltonian, {hj (t)} (Sec. B 2 a), and with
respect to the coefficients {cj } that specify the initial
K̇(t) = i [K(t), H(t)] (B4)
condition Eq. (16) (Sec. B 2 b).
=i aj (t)hk (t) [bj , bk ] (B5)
X

jk

= aj (t)hk (t)λklj bl . (B6)

X
a. With respect to system Hamiltonian
jkl

Comparing coefficients for bl (because of linear indepen- Let hk,m be constant within a time slice m, represent-
dence), gives ing a discretized time-dependence hk (t).
P In this time-
slice, the matrix G(t) is given by Gm = k hk,m Λk where
ȧl (t) = λklj aj (t)hk (t), (B7)
X
Λk is a d × d matrix built from structure coefficients λklj ,
jk
and d is the number of elements in the Lie algebra b.
which can be written compactly as With a total of M time-slices, where each time-slice du-
ration is τm = tm − tm−1 , the unitary propagator from
ȧ(t) = G(t)a(t), (B8) the solution of Eq. (B8), is given by Um = eGm τm . One
 15

can write, Let A = ∥aT ∥, V = ∥a(tf )∥ and F = aT · a(tf ), such that

J = 1 − F/AV . For each ℓ component in a(tf ),
aT · a(tf )
 
∂J ∂
=− (B12)
∂hk,m ∂hk,m ∥a(0)∥∥aT ∥ 1
 
∂J ∂ F
=− (B24)
a(tf ) ∂aℓ (tf ) A ∂aℓ (tf ) V
!
†
∂ aT
=− (B13)
1 V ∂ ℓ F − F ∂ℓ V
 
∂hk,m ∥a(0)∥∥aT ∥
=− , (B25)
†
! A V2
∂ aT UM UM −1
. . . U1 a(0)
=− (B14)
∂hk,m ∥a(0)∥∥aT ∥ where ∂ℓ = ∂/∂aℓ (tf ). With F = aT,j aj (tf ), one has
P
j
i 21
1 ∂  †
hP
2
∂ℓ F = aT,ℓ . Additionally with V = (t ) , one

=− aT UM UM −1 . . . U1 a(0) a
j j f
∥a(0)∥∥aT ∥ ∂hk,m
has
(B15)
1 ∂U − 21
=− (B16)
m
b†

∥a(0)∥∥aT ∥ m+1 ∂hk,m
am−1 1 X 2 aℓ (tf )
∂ℓ V = a (tf ) 2aℓ (tf ) = . (B26)
2 j j V
where b†m+1 = aT
†
UM UM −1 . . . Um+1 (backwards propa-
gation) and am−1 = Um−1 . . . U2 U1 a (0) (forwards prop-
Therefore, we have
agation).
For any differentiable matrix-valued function X(α),
∂J 1 V aT,ℓ − F [aℓ (tf ) /V ]
=− (B27)
∂aℓ (tf ) V2
Z 1
∂ X(α) ∂X sX(α) A
e = e(1−s)X(α) e ds. (B17)
∂α ∂α aT,ℓ F aℓ (tf )
0 =− + , (B28)
AV AV 3
For ∂Um /∂hk,m , one has Um = eGm τm and α = hk,m .
Thus ∂X/∂α becomes τm ∂Gm /∂hk,m = τm Λk , and one or in column vector form,
can write
Z 1 aT F
∂Um ∇a(tf ) J = − + a(tf ) = vT . (B29)
= τm e(1−s)Gm τm Λk esGm τm ds (B18) AV AV 3
∂hk,m 0
Let {i1 , . . . , in+1 } ⊂ {1, . . . , d} be the index set that be-
Z 1
= Um τm e−sGm τm Λk esGm τm ds (B19) longs to h, the maximal abelian set in m, where d is
0
Z τm the number of elements in the Lie algebra b. h contains
= Um e−tGm Λk etGm dt, (B20) n+1 elements. Let P be a selection matrix with elements
0 Pℓ,i = δℓ,hi , such that
where in the last line the substitution s = t/τn is made.
Taylor expanding the exponentials and computing the a(0) = P c, (B30)
integration gives,
where c = (c1 , . . . , cn+1 ) is a vector of coefficients for
⊤
τ2 the initial condition. Therefore, ∂a(0)/∂c = P ,
 
∂Um 3
= τm Λk Um − m [Gm , Λk ] Um + O τm


,
∂hk,m 2
(B21) ∂J ∂J ∂a(tf )
and finally, = (B31)
∂c ∂a(tf ) ∂c
∂J 1 τ2 =
∂J ∂a(0)
(B32)
≈− b†m+1 (τm Λk − m [Gm , Λk ])am ∂a(tf )
UM UM −1 . . . U1 .
∂hk,m ∥a(0)∥∥aT ∥ 2 ∂c
(B22)
Noting that ∂J /∂a(tf ) is the Jacobian row correspond-
ing to the column gradient ∇a(tf ) J , i.e., ∂J /∂a(tf ) =
b. With respect to initial condition (∇a(tf ) J )⊤ , we can write

The coefficients {ci } that specify the initial condition, ∂J ∂J ∂a(0)

= UM UM −1 . . . U1 (B33)
n ∂c ∂a(tf ) ∂c
K0 = cj Zj + cn+1 Z, (B23) = vT (B34)
X
⊤
UM UM −1 . . . U1 P
j=1
= v0⊤ P, (B35)
are additionally found during the optimization. Note
that cn+1 was given as c0 in the main text (Eq. (21)). where v0⊤ = vT
⊤
UM UM −1 . . . U1 .
 16

3. Further technicalities of the time-dependent factors {hj (t)}. At the cost of

longer computational running time, the optimization rou-
a. Objective function tine seems to escape unwanted local minima more effec-
tively. This is particularly useful for resolving the quan-
The operator infidelity J (K(tf ), KT ) (Eq. (17)) is min- tum speed limit, Figs. 5 and 10, in which optimal control
imal if and only if the operators K(tf ) and KT coincide solutions with operator infidelities of O(10−10 ) are typi-
up to a positive scalar factor, i.e. K(tf ) = cKT . To cally required. In this data, a discretization of N = 150n
achieve the actual goal, where K(tf ) and KT exactly co- was found to be sufficient, where optimal control solu-
incide for J (K(tf ), KT ) = 0, one can rescale the initial tions of sufficient accuracy to resolve the quantum speed
condition K0 . With, limit are often found with a single run of the algorithm.
This is the main reason why the data is presented only
f (c) = ∥K(tf ) − cKT ∥
2
(B36) up to moderate system sizes (up to n = 10), where the
2 2
computational running times at this discretization are
= ∥K(tf )∥ − 2c ⟨K (tf ) , KT ⟩ + c2 ∥KT ∥ , (B37) reasonable.
Secondly, the optimization routine can be run for a less
Given Hermitian operators, the expression is minimized optimal discretization, and then one can improve upon
by scalar the yielded solution. In particular, with a discretization
of N = 20n, the optimization will almost certainly con-
tr(K(tf )KT )
c= 2) , (B38) verge to an infidelity that is considerably higher than
tr(KT what is desired, especially for large system sizes n ≥ 20.
Taking this suboptimal solution as the initial guess for
satisfying f ′ (c) = 0. With this expression being effi-
the next optimization, with the same discretization, one
ciently evaluable using the coefficient vectors, Eq. (3) can
can iteratively improve on the original solution. Alter-
be restated as
natively, one can start with the optimal control solution
U (tf ) c−1 K0 U † (tf ) = c−1 K(tf ) ≈ KT (B39) for a nearby point (in the phase diagram), which should
have an easier optimization landscape, and then pertur-
enabling the required unitary mapping. batively progress towards the desired point. Namely, one
starts with the initial guess as the solution of the previ-
ous optimization and specifies the target as the perturbed
b. Discretization point.
In all cases, the only restriction to finding a better so-
lution is time. One can always improve upon the yielded
There are points in the phase diagram (Fig. 2) that
solution if time allows, as opposed to restrictions involv-
are more difficult to find desired optimal control solu-
ing memory (as is the case for explicit state-vector or
tions than others. These are typically the critical points,
density matrix control methods). To that end, a solu-
P1 , . . . , P4 . Because, regardless of the point in the phase
tion must always be accepted within a reasonable time
diagram, there is a guaranteed solution that maps the
frame. For the data presented in this paper, Figs. 4 and
initial thermal state ϱβ (0) to the target thermal state
9, operator infidelities of O(10−4 ) (resulting in worst-case
ϱβT via a unitary generated by the system Hamiltonian critical region thermal state infidelities between O(10−3 )
(Eq. (11)) (owing to Eq. (7)), the difficulty can be at- and O(10−2 )) are typically accepted solutions given the
tributed to a challenging optimization landscape rather accuracy of today’s quantum hardware.
than a fundamental barrier to reaching the solution.
To address such points, there are some helpful aspects
to consider. Firstly, one can select a finer discretization

[1] C. Jarzynski, Nonequilibrium equality for free energy dif- A. Gilyén, Quantum thermal state preparation (2023),
ferences, Phys. Rev. Lett. 78, 2690 (1997). arXiv:2303.18224 [quant-ph].
[2] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, [5] E. Y. Loh, J. E. Gubernatis, R. T. Scalettar, S. R. White,
From quantum chaos and eigenstate thermalization to D. J. Scalapino, and R. L. Sugar, Sign problem in the nu-
statistical mechanics and thermodynamics, Advances in merical simulation of many-electron systems, Phys. Rev.
Physics 65, 239 (2016). B 41, 9301 (1990).
[3] R. Sagastizabal, S. Premaratne, B. Klaver, M. Rol, [6] M. Troyer and U.-J. Wiese, Computational complexity
V. Negîrneac, M. Moreira, X. Zou, S. Johri, and fundamental limitations to fermionic quantum monte
N. Muthusubramanian, M. Beekman, et al., Variational carlo simulations, Phys. Rev. Lett. 94, 170201 (2005).
preparation of finite-temperature states on a quantum [7] R. Sweke, I. Sinayskiy, and F. Petruccione, Simulation
computer, npj Quantum Information 7, 130 (2021). of single-qubit open quantum systems, Phys. Rev. A 90,
[4] C.-F. Chen, M. J. Kastoryano, F. G. Brandão, and 022331 (2014).
 17

[8] B. Dive, F. Mintert, and D. Burgarth, Quantum simula- [25] A. Eckardt, Colloquium: Atomic quantum gases in pe-
tions of dissipative dynamics: Time dependence instead riodically driven optical lattices, Rev. Mod. Phys. 89,
of size, Phys. Rev. A 92, 032111 (2015). 011004 (2017).
[9] J. T. Barreiro, M. Müller, P. Schindler, D. Nigg, [26] P. Smacchia, L. Amico, P. Facchi, R. Fazio, G. Florio,
T. Monz, M. Chwalla, M. Hennrich, C. F. Roos, P. Zoller, S. Pascazio, and V. Vedral, Statistical mechanics of the
and R. Blatt, An open-system quantum simulator with cluster ising model, Phys. Rev. A 84, 022304 (2011).
trapped ions, Nature 470, 486 (2011). [27] R. Raussendorf, D. E. Browne, and H. J. Briegel,
[10] C.-F. Chen, M. J. Kastoryano, and A. Gilyén, An effi- Measurement-based quantum computation on cluster
cient and exact noncommutative quantum gibbs sampler states, Phys. Rev. A 68, 022312 (2003).
(2025), arXiv:2311.09207 [quant-ph]. [28] C. Ding, Phase transitions of a cluster ising model, Phys.
[11] P. Rall, C. Wang, and P. Wocjan, Thermal State Prepa- Rev. E 100, 042131 (2019).
ration via Rounding Promises, Quantum 7, 1132 (2023). [29] Eq. (3) might look deceptively similar to the evolution
[12] J. Guo, O. Hart, C.-F. Chen, A. J. Friedman, and A. Lu- A(t) = U (t)† A(0)U (t) of an observable A in the Heisen-
cas, Designing open quantum systems with known steady berg picture, all of the following discussion applies to the
states: Davies generators and beyond, Quantum 9, 1612 Schrödinger picture with time-independent observables.
(2025). [30] M. Orozco-Ruiz, N. H. Le, and F. Mintert, Quantum con-
[13] D. Hahn, S. Parameswaran, and B. Placke, Provably effi- trol without quantum states, PRX Quantum 5, 040346
cient quantum thermal state preparation via local driving (2024).
(2025), arXiv:2505.22816 [quant-ph]. [31] É. Cartan, Sur la structure des groups de transformations
[14] D. Hahn, R. Sweke, A. Deshpande, and O. Shtanko, Effi- finis et continus (Nony, 1894).
cient quantum gibbs sampling with local circuits (2025), [32] N. Khaneja and S. Glaser, Cartan decomposition of
arXiv:2506.04321 [quant-ph]. SU(2n ), constructive controllability of spin systems
[15] M. Hagan and N. Wiebe, The thermodynamic cost of ig- and universal quantum computing (2000), arXiv:quant-
norance: thermal state preparation with one ancilla qubit ph/0010100 [quant-ph].
(2025), arXiv:2502.03410 [quant-ph]. [33] E. Kökcü, T. Steckmann, Y. Wang, J. K. Freericks, E. F.
[16] J. Langbehn, G. Mouloudakis, E. King, R. Menu, Dumitrescu, and A. F. Kemper, Fixed depth hamiltonian
I. Gornyi, G. Morigi, Y. Gefen, and C. P. Koch, Uni- simulation via cartan decomposition, Phys. Rev. Lett.
versal cooling of quantum systems via randomized mea- 129, 070501 (2022).
surements (2025), arXiv:2506.11964 [quant-ph]. [34] S. Deffner and E. Lutz, Quantum speed limit for non-
[17] J. Lloyd and D. A. Abanin, Quantum thermal state markovian dynamics, Phys. Rev. Lett. 111, 010402
preparation for near-term quantum processors (2025), (2013).
arXiv:2506.21318 [quant-ph]. [35] S. Deffner and S. Campbell, Quantum speed limits: from
[18] M. Scandi and Á. M. Alhambra, Thermalization in open heisenberg’s uncertainty principle to optimal quantum
many-body systems and kms detailed balance (2025), control, Journal of Physics A: Mathematical and The-
arXiv:2505.20064 [quant-ph]. oretical 50, 453001 (2017).
[19] O. Shtanko and R. Movassagh, Preparing thermal states [36] W. Cottrell, B. Freivogel, D. M. Hofman, and S. F.
on noiseless and noisy programmable quantum processors Lokhande, How to build the thermofield double state,
(2023), arXiv:2112.14688 [quant-ph]. Journal of High Energy Physics 2019, 1 (2019).
[20] Z. Ding, Y. Zhan, J. Preskill, and L. Lin, End-to-end [37] J. Wu and T. H. Hsieh, Variational thermal quantum
efficient quantum thermal and ground state preparation simulation via thermofield double states, Phys. Rev. Lett.
made simple (2025), arXiv:2508.05703 [quant-ph]. 123, 220502 (2019).
[21] K. Temme, T. J. Osborne, K. G. Vollbrecht, D. Poulin, [38] R. Wiersema, E. Kökcü, A. F. Kemper, and B. N.
and F. Verstraete, Quantum metropolis sampling, Nature Bakalov, Classification of dynamical lie algebras of 2-
471, 87 (2011). local spin systems on linear, circular and fully connected
[22] M. J. Kastoryano and K. Temme, Quantum logarithmic topologies, npj Quantum Information 10, 110 (2024).
sobolev inequalities and rapid mixing, Journal of Math- [39] L. Stefanescu, L. Edwards-Pratt, J. O’Connor,
ematical Physics 54 (2013). E. Tsegaye, N. H. Le, and F. Mintert, Robust im-
[23] S. P. Jordan and E. Farhi, Perturbative gadgets at arbi- plicit quantum control of interacting spin chains, Phys.
trary orders, Phys. Rev. A 77, 062329 (2008). Rev. A 112, 012609 (2025).
[24] I. M. Georgescu, S. Ashhab, and F. Nori, Quantum sim- [40] A. van Lomwel, Data for van lomwel et al., fast
ulation, Rev. Mod. Phys. 86, 153 (2014). thermal state preparation beyond native interactions.,
10.5281/zenodo.17936714 (2026).