# Deep reinforcement learning for preparation of thermal and prethermal
quantum states

Shotaro Z. Baba<EMAIL_ADDRESS>Department of Applied Physics, The
University of Tokyo, Hongo, Bunkyo-ku, Tokyo 113-8656, Japan  Nobuyuki
Yoshioka<EMAIL_ADDRESS>Department of Applied Physics, The
University of Tokyo, Hongo, Bunkyo-ku, Tokyo 113-8656, Japan  Theoretical
Quantum Physics Laboratory, RIKEN Cluster for Pioneering Research (CPR), Wako-
shi, Saitama 351-0198, Japan  JST, PRESTO, 4-1-8 Honcho, Kawaguchi, Saitama,
332-0012, Japan  Yuto Ashida  Department of Physics, The University of Tokyo,
Hongo, Bunkyo-Ku, Tokyo 113-0033, Japan  Takahiro Sagawa  Department of
Applied Physics, The University of Tokyo, Hongo, Bunkyo-ku, Tokyo 113-8656,
Japan  Quantum-Phase Electronics Center (QPEC), The University of Tokyo, Tokyo
113-8656, Japan

###### Abstract

We propose a method based on deep reinforcement learning that efficiently
prepares a quantum many-body pure state in thermal or prethermal equilibrium.
The main physical intuition underlying the method is that the information on
the equilibrium states can be efficiently encoded/extracted by focusing only
on a few local observables, relying on the typicality of equilibrium states.
Instead of resorting to the expensive preparation protocol that adopts global
features such as the quantum state fidelity, we show that the equilibrium
states can be efficiently prepared only by learning the expectation values of
local observables. We demonstrate our method by preparing two illustrative
examples: Gibbs ensembles in non-integrable systems and generalized Gibbs
ensembles in integrable systems. Pure states prepared solely from local
observables are numerically shown to successfully encode the macroscopic
properties of the equilibrium states. Furthermore, we find that the
preparation errors, with respect to the system size, decay exponentially for
Gibbs ensembles and polynomially for generalized Gibbs ensembles, which are in
agreement with the finite-size fluctuation within thermodynamic ensembles. Our
method paves a new path to studying thermodynamic and statistical properties
of quantum many-body systems in quantum hardware.

##  I Introduction

Preparation of a desired quantum many-body state is an essential task that
plays a significant role in quantum computing [1, 2], quantum metrology [3],
and quantum communication [4]. Specifically, the thermal state is one of the
most important targets in quantum state preparation tasks [5, 6, 7, 8, 9, 10,
11] from both theoretical and experimental viewpoints. A common strategy is to
employ numerical methods such as CRAB [12, 13], GRAPE [14], and Krotov [15].
However, it must be noted that all these methods suffer from the exponential
growth of computational cost, and furthermore require detailed knowledge about
nonequilibrium properties of the system. Therefore, it is desirable to
construct a preparation protocol that employs only a little knowledge during
the learning task.

A surging technology to extract the essential feature in quantum systems with
prohibitively large exploration space is machine learning, which has
exemplified its capacity in a wide range of physics [16, 17, 18, 19, 20, 21,
22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37]. Successful
applications include representations of quantum many-body states with neural
networks (NNs) [18, 19, 38, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 39],
quantum state tomography [31, 32, 33], and phase classification [34, 35, 40],
to name a few. In particular, a branch of machine learning called
reinforcement learning (RL) [41] has been recognized as a powerful tool to
perform quantum state preparation [42, 43, 44, 45, 46, 47, 48, 49, 50]. The RL
is designed to discover an efficient policy that maximizes a given reward
through trial-and-error learning on the behavior of the environment. Several
previous studies utilize the algorithm that adopts the deep RL framework; the
quantitative evaluation of the action, or the reward, determined by the
algorithm makes full use of the capability of NNs to approximate the high-
dimensional nonlinear functions [50, 51, 48, 47, 46, 45, 49, 52]. While a bulk
of previous works choose fidelity as the reward, its computation for quantum
many-body systems requires exponentially large resources in either numerics or
experiment, and thus fidelity is not practical to scale up.

In this work, we propose a deep-RL-based method that only relies on local
measurements to prepare thermal and prethermal pure states described by Gibbs
and generalized Gibbs ensembles (GGEs) [53, 54, 55, 56]. The underlying
physical intuition is that, we may take advantage of the typicality of
equilibrium states [57, 58, 7, 59] to prepare them using solely local
observables, and not rely on global features such as fidelity. We numerically
find that, although the deep RL agent is only informed of the local
information on the thermodynamic ensembles, the accuracy of the prepared state
improves exponentially with the system size for Gibbs ensembles, whereas the
improvement is polynomial for GGEs.

The remainder of the paper is organized as follows. In Sec. II, we present an
overview of the framework of the deep reinforcement learning. The application
of reinforcement learning to quantum state preparation is described in Sec.
III, which includes the core proposal of our work, i.e., the local preparation
of thermal and prethermal pure quantum states leveraging the typicality of
equilibrium states. After the framework is presented, we give the numerical
demonstration for preparation of equilibrium states described by Gibbs
ensembles and GGEs in Sec. IV and V, respectively. Finally, we give the
conclusion and discussion in Sec. VI.

##  II Reinforcement learning

The general framework of the reinforcement learning (RL) can be concisely
expressed as a procedure to train an agent how to interact with the
environment through optimization of cumulative reward [see Fig. 1]. Namely,
the goal is to discover the best policy $\pi$ that outputs the sequence of
actions $\mathcal{A}=\\{a_{t}\\}_{t}$ based on observations of the environment
$\mathcal{O}=\\{o_{t}\\}_{t}$, such that the feedback realizes the most
desired behavior quantified by the rewards $\mathcal{R}=\\{r_{t}\\}_{t}$.
Typically, all the events are discretized so that each value can be well-
defined at each time step $t$.

A practical strategy widely used in the community of machine learning is
Q-learning. Here, one aims to find the best approximation of the optimal
action-value function as [41]

$\displaystyle
Q^{*}(o,a)=\max_{\pi}\mathbb{E}_{\pi}\left[r_{t}+\sum_{n=1}^{\infty}\gamma^{n}r_{t+n}\middle|o_{t}=o,\
a_{t}=a,\pi\right],$ (1)

which is the maximum sum of rewards $r_{t}$ discounted by $\gamma$
($0<\gamma<1$) in a stochastic policy that chooses action according to some
probability distribution as $\pi\left\lparen a|o\right\rparen=\Pr\left\lparen
a|o\right\rparen$. A powerful flavor of Q-learning uses the deep NNs to
represent the action-value function, and hence referred to as the deep RL
algorithms [60, 61]. The extraordinary representative power of NNs have been
found to achieve successful applications of the deep RL algorithms in numerous
fields that are not necessarily limited to computer science but also natural
science, materials sciences, and so on.

In this paper, we focus on a non-distributed implementation [62] of a deep RL
algorithm called R2D2 [63]. R2D2 is a type of deep Q-learning algorithm [64],
and assumes the agent can obtain partial information about the state of the
environment. As we show the architecture in Appendix A, the NN used in R2D2
includes a Long Short-Term Memory (LSTM) layer, and therefore the action-value
function $Q$ computed by the NN at step $t$ depends not only on the instant
observation $o_{t}$, but also on the previous observations
${\left\\{o_{t^{\prime}}\right\\}}_{t^{\prime}\leq t}$ [65]. This feature
enables the NN to handle time-series inputs and develop the capability in a
partially observed environment.

##  III Reinforcement learning for quantum state preparation

###  III.1 Global state preparation

Next, we review the general protocol to prepare a desired isolated quantum
many-body state using the deep RL framework. Concretely, we aim to prepare a
target quantum state $\ket{\psi_{\rm target}}$ from an easily prepared quantum
state $\ket{\psi_{0}}$, assuming that a set of unitary $\\{U_{i}\\}_{i}$ is
available at any time step. By finding the best sequence of unitaries, we try
to approximate the target state as

$\ket{\psi_{\rm target}}\approx\prod_{t}U_{i_{t}}\ket{\psi_{0}}.$ (2)

It is straightforward to see that such a problem setup is in a great
connection with the RL; we identify the quantum many-body system with the
environment and the available set of unitaries $\\{U_{i}\\}_{i}$ with the
action candidates $\\{a_{t}\\}$ at each time step.

Regarding the observation $o_{t}$, many works have proposed to use the results
for the measurements on the target system [45, 46, 47]. Meanwhile, when the
both initial and target states are fixed during the whole training, we can
expect that the action history $\\{a_{t^{\prime}}\\}_{t^{\prime}\leq t}$
contains enough information to find out the desired protocol [43, 44].

As for the reward $r_{t}$, numerous existing works have considered global
features such as the fidelity
$F(\rho_{t},\rho_{\mathrm{target}})={\left\lparen\operatorname{Tr}\sqrt{\sqrt{\rho_{t}}\rho_{\mathrm{target}}\sqrt{\rho_{t}}}\right\rparen}^{2}$
[42, 43, 48, 50], where $\rho_{t}$ is the density operator of the controlled
system at time step $t$. We hereafter refer to such protocols as global state
preparation protocols. These methods have successfully prepared ground states
of quantum many-body spin systems [42], metastable states of the quantum
Kapitza oscillator [43], and highest excited states of multi-level quantum
systems [50].

Figure 1:  A graphical illustration of the deep RL framework to prepare
thermal and prethermal pure quantum states for thermodynamic ensembles. The
agent represented by the deep NN is trained to discover a policy that
maximizes the cumulative reward, which is computed based on the local
observables of the quantum many-body system (environment). At each time step,
the NN takes the result of observation as input and outputs the action-value
function, from which the action (or the control operation) on the environment
is determined. In this work, the action is chosen from a set of unitaries that
are presumably available in the preparation protocol.

Having bridged between the notions in quantum control and the RL, we can train
the NN-based agent to find the best control on the quantum system via
searching for approximation of the optimal action-value function $Q^{*}(o,a)$
(1). Note that the evolution of the quantum state may either be experimentally
implemented or numerically simulated, as long as the reward function for the
deep RL agent can be readily obtained. Once the training is completed, we
determine the preparation protocol as $\\{a^{*}_{t}\\}_{t}$ by choosing the
actions so that $a^{*}_{t}=\operatorname*{arg\,max}_{a}Q(o_{t},a)$ at each
time step $t$, which shall maximize the reward.

###  III.2 Local state preparation

Now we are ready to describe the preparation protocol for thermal and
prethermal pure quantum states that solely relies on local observables,
instead of querying for costly global features such as the fidelity. In the
following, as opposed to the global state preparation, we refer to the
following scheme as local state preparation.

The central idea of the local state preparation is to leverage the typicality
of pure quantum states [57, 58, 7, 59]. The typicality refers to the fact that
the overwhelming majority of pure quantum states with the same local conserved
quantities are indistinguishable by local observables. Such states are also
deemed to be macroscopically indistinguishable. It is natural to expect that,
by utilizing the typicality, we can prepare a state that encodes the
macroscopic property of the equilibrium states by solely controlling local
observables. Specifically, we perform local preparation by making the
expectation values of local observables close to the equilibrium state, and
then letting the system relax to equilibrium through a unitary evolution
without any control.

We may take various forms of reward to learn the local observables of the
target ensemble. In this work, we formulate the reward function as the inverse
of the deviations of the expectation values of local observables from that of
the target states:

$\displaystyle r_{t}\coloneqq\frac{1}{\left\lvert
M_{t}-M_{\mathrm{target}}\right\rvert+\varepsilon},$ (3)

where $M_{t}=(\braket{O_{1}}_{t},\braket{O_{2}}_{t},\ldots)$ and
$M_{\mathrm{target}}=(\braket{O_{1}}_{\mathrm{target}},\braket{O_{2}}_{\mathrm{target}},\ldots)$
are vectors consisting of expectation values $\braket{O}_{t({\rm
target})}:={\rm Tr}[\rho_{t({\rm target})}O]$. The small constant
$\varepsilon$ is also introduced to prevent divergence of the reward function.

As a concrete target for the demonstration of our local preparation protocol,
we choose Gibbs ensembles and GGEs as illustrative examples. The evolution of
the quantum many-body state is numerically simulated in an exact manner, while
in principle we may also employ approximate methods that rely on, e.g., a
variational representation such as tensor networks or neural networks. One may
alternatively implement the proposed protocol directly on experimental device
as well. In the following, we proceed to describe the detailed properties of
the thermodynamic ensembles and the expected preparation efficiency under the
presence of typicality.

####  III.2.1 Gibbs ensemble for non-integrable systems

Let us consider a non-integrable system with the energy being the only local
conserved quantity such that its equilibrium is described by a Gibbs ensemble.
Relatedly, pure states belonging to a given microcanonical shell of the system
share their macroscopic properties. This class of typicality is referred to as
canonical typicality. One of the most prominent example of canonical
typicality can be illustrated under the Haar measure on the space of pure
states under some constraint $R$ (e.g., energy) as [58]

$\displaystyle\braket{\left\lVert\rho_{A}-\Omega_{A}\right\rVert_{1}}\leq\sqrt{\frac{d_{A}^{2}}{d_{R}}},$
(4)

where $\rho_{A}=\operatorname{Tr}_{\bar{A}}\left[\ket{\psi}\bra{\psi}\right]$
is the reduced density operator obtained by tracing out the the complement of
subsystem $A$ for $\ket{\psi}\bra{\psi}$ with Hilbert space dimension $d_{A}$.
On the other hand,
$\Omega_{A}=\operatorname{Tr}_{\bar{A}}\left[\mathds{1}_{R}/d_{R}\right]$ is
the reduced density operator for the projection operator $\mathds{1}_{R}$,
i.e., the maximally mixed state in the Hilbert subspace under the constraint
$R$ whose dimension is $d_{R}$. Note that the bracket $\braket{\cdot}$
concerns the average regarding the Haar measure on the constrained Hilbert
space, and $\left\lVert
A\right\rVert_{1}=\operatorname{Tr}\sqrt{A^{\dagger}A}$.

Equation (4) means that the average distance between a randomly chosen pure
state and the maximally mixed state in the constrained Hilbert space decays
polynomially with the Hilbert space size as $d_{R}^{-0.5}$, that is,
exponentially with the system size in general quantum systems. The
indistinguishability of the pure quantum state from the microcanonical
ensemble leads us to expect that, once a state is prepared to be within the
target energy shell, the prepared state captures the macroscopic properties
with an accuracy that improves exponentially with the system size 111 Gibbs
ensembles, which is our target, differs from microcanonical ensembles.
Nevertheless, a microcanonical ensemble is locally equivalent to a Gibbs
(canonical) ensemble for translation invariant short-ranged Hamiltonians, even
if the system is finite [82, 83]. In other words, the microcanonical and the
canonical expectation values almost coincide when we look at the subsystem
that is not too large..

We emphasize that the key of local preparation protocol is to encode the
prepared state into the target energy shell, which requires more than merely
learning the expectation value of the energy. In other words, even if the
expectation value is correctly learned, the prepared state may correspond to a
superposition of pure states that belong to other energy shells. Such a
situation may cause deviation in other physical observables. In this work, we
attempt to address this problem by letting the RL agent learn other
macroscopic observables in addition to energy. It is in fact highly nontrivial
to determine how many observables we need to embed the prepared state into the
energy shell. We find that, for the non-integrable transverse-field Ising
chain, it suffices to take only the total magnetization
$\sum_{l}\sigma_{l}^{z}/L$ [see Sec. IV.2 for the numerical demonstration].

As another possibly effective method to assure the prepared state to be
encoded in the desired energy shell, we propose to incorporate the variance of
observables $\braket{O^{2}}_{\rm target}-\braket{O}_{\rm target}^{2}$ into the
reward. For instance, it is obvious that the energy variance is suppressed
when the prepared state is in the target energy shell. This is actually not
limited to the energy; thanks to the typicality, we may employ any macroscopic
observable for this purpose.

####  III.2.2 Generalized Gibbs ensemble for integrable systems

As opposed to the non-integrable systems, integrable systems have an extensive
number of conserved quantities, which are also called integrals of motion
(IOM) in the literature. The equilibrium states in integrable systems are
known to be described rather by GGEs [54]:

$\displaystyle\rho_{\mathrm{GGE}}=\frac{\exp\left[-\sum_{m}\lambda_{m}\hat{\mathcal{I}}_{m}\right]}{\operatorname{Tr}\left[\exp\left[-\sum_{m}\lambda_{m}\hat{\mathcal{I}}_{m}\right]\right]},$
(5)

where $\\{\hat{\mathcal{I}}_{m}\\}$ is the full set of the IOMs, and
$\left\\{\lambda_{m}\right\\}$ is the corresponding set of the Lagrange
multipliers which dictates the distribution over expectation values of the
IOMs.

To discuss the typicality in the set of pure states with close expectation
values of IOM, the authors of Ref. [67] introduced a notion of a statistical
ensemble named the generalized microcanonical ensemble (GME). In parallel to
the ordinary microcanonical ensemble, the GME is constructed by assigning
equal weight to all eigenstates whose IOMs are close to some certain values
that identify the ensemble. It has been pointed out Ref. [67] that the
standard deviations of local observables within such “a window of IOM” decay
polynomially as

$\displaystyle\sigma_{\mathrm{\rm loc}}\propto L^{-0.5},$ (6)

where $L$ is the system size. Therefore, by following a parallel discussion as
in the case for Gibbs ensembles in non-integrable systems, we may also expect
that the local preparation protocol works for GGEs in integrable systems as
well, with its accuracy improving polynomially with the system size.

Let us remark on another supporting argument based on the truncation of GGEs
itself. While Eq. (5) takes all possible IOMs into account, we expect that the
macroscopic behavior in terms of local observables can be extracted by
considering local conserved quantities. As such, here, we aim to capture the
truncated alternative of the statistical ensemble by focusing on the local
integrals of motion (LIOM). Concretely, we denote the LIOMs that acts at most
$n+1$ neighboring sites as $\hat{I}_{n}^{\sigma}$ with $\sigma$ denoting some
additional label, and introduce a locality-constrained variant of GGE which is
known as the truncated GGE (tGGE) [68]:

$\displaystyle\rho_{\mathrm{tGGE},n_{\mathrm{local}}}=\frac{\exp\left[-\sum_{n=0}^{n_{\mathrm{local}}}\sum_{\sigma}\lambda_{n}^{\sigma}\hat{I}_{n}^{\sigma}\right]}{\operatorname{Tr}\left[\exp\left[-\sum_{n=0}^{n_{\mathrm{local}}}\sum_{\sigma}\lambda_{n}^{\sigma}\hat{I}_{n}^{\sigma}\right]\right]},$
(7)

which only includes the LIOMs with $n\leq n_{\mathrm{local}}$. It is natural
to expect that a tGGE gives a good approximation of $\rho_{\rm GGE}$ in terms
of local quantities, and furthermore
$\operatorname{Tr}_{\bar{A}}\left[\rho_{\mathrm{tGGE},n_{\mathrm{local}}}\right]\approx\operatorname{Tr}_{\bar{A}}\left[\rho_{\mathrm{GGE}}\right]$.
For instance, Ref. [68] has investigated the transverse field Ising chain in
the integrable regime, and found that tGGEs approximates the corresponding
GGEs when $n_{\rm local}$ is larger than the size of subsystem $A$.

We remark that the local preparation protocol for the integrable systems aims
to construct a tGGE rather than the original GGE. In this sense, we expect
that the validity is not assured for observables with higher $n_{\rm local}$,
for which the discrepancy between the GGE and tGGE is non-negligible. We
discuss this point more in detail in Sec. V.2.

##  IV Application to Gibbs ensembles

###  IV.1 Model and setup

As a demonstration for local preparation of thermal pure states described by
Gibbs ensembles, we consider the transverse field Ising model on a chain with
the periodic boundary condition:

$\displaystyle\hat{H}_{\mathrm{Ising}}=\sum_{l=1}^{L}\left[J\hat{\sigma}_{l}^{z}\hat{\sigma}_{l+1}^{z}+h\hat{\sigma}_{l}^{z}+g\hat{\sigma}_{l}^{x}\right],$
(8)

where $\hat{\sigma}_{l}^{x},\ \hat{\sigma}_{l}^{y},\ \hat{\sigma}_{l}^{z}$ are
the Pauli operators at site $l$, $L$ is the system size, $J$ is the amplitude
of the Ising interaction, and $h~{}(g)$ is the strength of the longitudinal
(transverse) magnetic field. In the following, the parameters are fixed as
$J=1,\ h=0.8090,\ g=0.9045$ so that the system is non-integrable [69]. As a
target state, we aim to prepare a thermal pure quantum state corresponding to
the Gibbs ensemble with inverse temperature $\beta=0.2$, where the initial
state before any quantum control is taken to be a product state
$\Ket{\downarrow\downarrow\cdots\downarrow}$. The total preparation time is
fixed as $T=24$ with the time step set as $\delta t=0.1$, which also
determines the total time step to be $240$.

The set of action candidates $\\{a_{t}\\}$ available for the RL agent is given
as $\\{e^{-iG_{k}\delta t}\\}_{k=1}^{6}$, where the time evolution generator
$G_{k}$ is chosen from the following six operators:

$\displaystyle\begin{split}&\hat{H}_{\mathrm{Ising}},\
\sum_{l=1}^{L}J\hat{\sigma}_{l}^{z}\hat{\sigma}_{l+1}^{z}+h\hat{\sigma}_{l}^{z},\
g\sum_{l=1}^{L}\hat{\sigma}_{l}^{x},\ \sum_{l=1}^{L}\hat{\sigma}_{l}^{y},\\\
&\sum_{l=1}^{L}\hat{\sigma}_{l}^{x}\hat{\sigma}_{l+1}^{y}+\hat{\sigma}_{l}^{y}\hat{\sigma}_{l+1}^{x},\
\sum_{l=1}^{L}\hat{\sigma}_{l}^{y}\hat{\sigma}_{l+1}^{z}+\hat{\sigma}_{l}^{z}\hat{\sigma}_{l+1}^{y}.\end{split}$
(9)

These terms are also used in Ref. [44], which discusses how to accelerate the
preparation of the ground state using counter-diabatic driving. Regarding the
local reward $r_{t}$, we include the energy density
$L^{-1}\hat{H}_{\mathrm{Ising}}$ and the magnetization density
$L^{-1}\sum_{l}\hat{\sigma}_{l}^{z}$. Note that these are the sum of local
operators acting on neighboring two sites at most.

###  IV.2 Numerical results

Figure 2:  (a) The learning curves for local preparation of the Gibbs ensemble
of transverse field Ising model in the non-integrable regime. The blue and red
curves denote the total reward and the deviation of the energy density
computed for the prepared states, respectively. The horizontal axis shows the
number of parameter updates. (b)–(e) The dynamics of local observables
generated by the state preparation protocol learned by the deep RL, with the
inset showing the absolute error from the corresponding expectation values of
the target Gibbs state (red line). Each panel displays (b) energy density
$\hat{H}_{\mathrm{Ising}}/L$, (c)
$Z^{(1)}/L=L^{-1}\sum_{l}\hat{\sigma}_{l}^{z}$, (d)
$Z^{(2)}/L=L^{-1}\sum_{l}\hat{\sigma}_{l}^{z}\hat{\sigma}_{l+1}^{z}$, and (e)
$Z^{(3)}/L=L^{-1}\sum_{l}\hat{\sigma}_{l}^{z}\hat{\sigma}_{l+1}^{z}\hat{\sigma}_{l+2}^{z}$.
Note that local observables in (b) and (c) are both included in the reward for
the deep RL agent. We observe that, not only the observables which are
included in the RL reward converge to the target value as in (b, c), but also
local observables absent in the reward readily converge to the target values
as in (d, e). For all plots, the system size is $L=16$ with the inverse
temperature set as $\beta=0.2$.

Now we proceed to present numerical results that successfully prepare thermal
pure quantum states, which encode both local observables used during the
training and also more non-local ones that are excluded from the reward
function. Figure 2 (a) shows an example of the learning curve of the RL agent
for $L=16$. We can see that, as the number of training episodes increases, the
RL agent learns the better protocols that achieve higher total reward and
smaller energy deviation. [See Appendix. C for the concrete protocol.]

Figures 2 (b)–(e) display the dynamics of local observables obtained by the
preparation protocol learned by the RL agent. Figures 2 (b) and (c) show the
observables $\hat{H}_{\mathrm{Ising}}$ and $\sum_{l}\hat{\sigma}_{l}^{z}$,
respectively, which are used for the reward function. Both of them converge to
the corresponding values of the Gibbs ensemble represented by the red
horizontal lines. It is more intriguing to see the convergent behaviors in
Figs. 2 (d) and (e), which show the dynamics of local observables that are not
included in the reward function, namely the two-point correlator
$\sum_{l}\hat{\sigma}_{l}^{z}\hat{\sigma}_{l+1}^{z}$ and three-point
correlator
$\sum_{l}\hat{\sigma}_{l}^{z}\hat{\sigma}_{l+1}^{z}\hat{\sigma}_{l+2}^{z}$,
respectively.

Figure 3:  Distance of reduced density operator between the target Gibbs state
and the prepared state, which are averaged over time and random training
instances. (a) Error of distance function $\bar{D}$ at various subsystem sizes
$LA$. The blue circles, orange downward triangles, green crosses, and red
upward triangles denote the data for $L=10,12,14,16$, respectively. (b) Decay
of error in the reduced density operator along with the corresponding energy
shell dimension $d$, i.e., the number of the eigenstates in the energy shell.
The dotted lines are guides to the eye showing the fits with $D=ad^{-b}$,
where powers $b$ are summarized in Table 3 of Appendix. B. The error bars
correspond to the standard deviations of different protocols learned
independently with random seeds. In both plots, the system is controlled by
the deep RL agent at $0\leq t\leq 24$, and then undergoes free evolution
$\exp[-i\hat{H}_{\mathrm{Ising}}t]$ in $24<t\leq 48$. Note that the distance
$\bar{D}$ concerns average over the free evolution period and also the random
training instance.

It is natural to wonder how accurately the general observables, or the reduced
density operator of local systems, are captured by the local preparation
protocol. For this purpose, we evaluate the distance between reduced density
operators of the prepared pure state and the target Gibbs ensemble as

$\displaystyle
D(\rho,\rho^{\prime})\coloneqq\frac{\left\lVert\rho-\rho^{\prime}\right\rVert_{\mathrm{F}}}{\sqrt{\left\lVert\rho\right\rVert_{\mathrm{F}}^{2}-\left\lVert\rho^{\prime}\right\rVert_{\mathrm{F}}^{2}}},$
(10)

where $\left\lVert
A\right\rVert_{\mathrm{F}}\coloneqq\sqrt{\operatorname{Tr}\left[A^{\dagger}A\right]}$
is the Frobenius norm 222Note that the chosen distance function can be
efficiently computed for Slater determinants, which is used to express free-
fermionic states obtained from Jordan-Wigner transformation. See Ref. [68] for
the detailed property of this distance. .

We can see from Fig. 3 (a) that the time average of the distance function
(10), which we denote as $\bar{D}$, is suppressed for the smaller subsystem
size $LA$. What is more interesting is the scaling with respect to the total
system size; we observe nontrivial suppression of error along with the
increase of the system size $L$.

We further zoom into the scaling of the error suppression, by performing
finite-size scaling on the averaged distance function as shown in Fig. 3 (b).
Here, we observe that the scaling of the suppression is given as

$\bar{D}=\mathcal{O}(d^{-b}),$ (11)

where $d$ is the dimension of the corresponding energy shell. To be concrete,
$d$ is obtained by counting the number of eigenstates included within the
energy window
$[\braket{\hat{H}_{\mathrm{Ising}}}_{\beta}-\varepsilon_{\mathrm{e}}L,~{}\braket{\hat{H}_{\mathrm{Ising}}}_{\beta}]$,
where $\braket{\hat{H}_{\mathrm{Ising}}}_{\beta}$ is the energy expectation
value of the target Gibbs ensemble at inverse temperature $\beta$, and the
shell width is fixed as $\varepsilon_{\mathrm{e}}=0.5$ [see Appendix D for
further discussion on the choice of $\varepsilon_{\mathrm{e}}$].

The scaling of the distance (11) implies that $\bar{D}$ decays exponentially
with the system size $L$. We argue that this is compatible with the scaling of
canonical typicality (4). This feature supports an expectation that our local
preparation protocol for the Gibbs ensembles becomes exponentially more
precise as the system size increases.

As a technical remark, we mention that the dimension shown in the figure
corresponds to the size of the symmetry-resolved Hilbert space, namely the
parity and the momentum.

##  V Application to generalized Gibbs ensemble

###  V.1 Model and setup

We next describe an even more intriguing system that is integrable so that the
thermodynamic ensembles is described by GGEs. Concretely, we consider the XX
model on a periodic chain with a longitudinal field:

$\displaystyle\hat{H}_{\mathrm{XX}}=-\sum_{l=1}^{L}\left[\frac{J}{2}\left(\hat{\sigma}_{l}^{x}\hat{\sigma}_{l+1}^{x}+\hat{\sigma}_{l}^{y}\hat{\sigma}_{l+1}^{y}\right)+h\hat{\sigma}_{l}^{z}\right],$
(12)

where $J$ is the amplitude of interaction, $h$ is the strength of the magnetic
field along the $z$-axis, and $L$ is the system size. In the following, the
parameters are fixed as $J=1$ and $h=2$.

The integrability of the XX model can be verified straightforwardly by mapping
into a non-interacting fermionic system via the Jordan-Wigner transformation:

$\displaystyle\hat{H}_{\mathrm{tb}}$ $\displaystyle=$
$\displaystyle-\sum_{l=1}^{L}[J(\hat{a}_{l}^{\dagger}\hat{a}_{l+1}+\hat{a}_{l+1}^{\dagger}\hat{a}_{l})+2h(\hat{a}_{l}^{\dagger}\hat{a}_{l}-\frac{1}{2})],$
(13) $\displaystyle=$ $\displaystyle-\sum_{k}\left[2\left\lparen J\cos
k+h\right\rparen\hat{n}_{k}-h\right],$ (14)

where $\hat{a}_{l}^{\dagger}$ and $\hat{a}_{l}$ are fermionic creation and
annihilation operators that are related to the spin operators as
$\hat{a}_{l}^{({\dagger})}=\prod_{l^{\prime}<l}\hat{\sigma}_{l^{\prime}}^{z}\hat{\sigma}_{l}^{+(-)}$
where $\hat{\sigma}_{l}^{\pm}=(\hat{\sigma}_{l}^{x}\pm
i\hat{\sigma}_{l}^{y})/2$. In the second row (14), we moved into the Fourier
space by introducing the mode occupation operator
$\hat{n}_{k}=\hat{a}_{k}^{\dagger}\hat{a}_{k}$ for
$\hat{a}_{k}=L^{-1/2}\sum_{l}\hat{a}_{l}e^{-ikl}$. We can construct LIOMs as
many as the system size $L$ by taking the linear combination as

$\displaystyle\hat{I}_{n}^{+}$
$\displaystyle=-2J\sum_{k}\cos{(nk)}\hat{a}_{k}^{\dagger}\hat{a}_{k}=-J\sum_{l}\left\lparen\hat{a}_{l}^{\dagger}\hat{a}_{l+n}+\hat{a}_{l+n}^{\dagger}\hat{a}_{l}\right\rparen,$
$\displaystyle\hat{I}_{n}^{-}$ $\displaystyle=-2J\sum_{k}\sin{\left\lparen
nk\right\rparen}\hat{a}_{k}^{\dagger}\hat{a}_{k}=iJ\sum_{l}\left\lparen\hat{a}_{l}^{\dagger}\hat{a}_{l+n}-\hat{a}_{l+n}^{\dagger}\hat{a}_{l}\right\rparen.$

For the convenience in the later discussion, we mention the rightmost sides to
remark that $\hat{I}_{n}^{\pm}$ can be explicitly expressed as sums over
hopping terms between the $n$-th nearest neighboring sites in the fermionic
picture.

As the target thermodynamic ensemble, we aim for the GGE constructed from the
LIOMs as

$\displaystyle\rho_{\mathrm{GGE}}=\frac{\exp\left[-\sum_{n}\sum_{\sigma=\pm}\lambda_{n}^{\sigma}\hat{I}_{n}^{\sigma}\right]}{\operatorname{Tr}\left[\exp\left[-\sum_{n}\sum_{\sigma=\pm}\lambda_{n}^{\sigma}\hat{I}_{n}^{\sigma}\right]\right]},$
(15)

where $\\{\lambda_{n}^{\sigma}\\}$ are the Lagrange multipliers. More
precisely, we set the initial state of the control system to be the ground
state of $\hat{H}_{\mathrm{tb}}$ in the space of the total particle number
$N_{\mathrm{f}}$, and aim to prepare a prethermal pure state corresponding to
a given set of Lagrange multipliers $\\{\lambda_{n}^{\sigma}\\}$ [see Appendix
E for details regarding the parameter choice of $\\{\lambda_{n}^{\sigma}\\}$].
The total preparation time is fixed as $T=40$ and the time step as $\delta
t=0.2$, thus the total number of time steps is $200$.

In parallel with the case for Gibbs ensembles, we allow the deep RL agent to
choose $a_{t}$ as an appropriate unitary from a set $\\{e^{-iG_{k}\delta
t}\\}_{k=1}^{15}$, where the time evolution generator $G_{k}$ is chosen from
the following fifteen operators:

$\displaystyle\begin{split}&\hat{H}_{\mathrm{XX}},~{}\sum_{l=1}^{L}\hat{a}_{l}^{\dagger}\hat{a}_{l}\cos\frac{m\pi}{L}l,~{}\sum_{l=1}^{L}\hat{a}_{l}^{\dagger}\hat{a}_{l}\sin\frac{m\pi}{L}l,\\\
&\sum_{l=1}^{L}\left\lparen\hat{a}_{l}^{\dagger}\hat{a}_{l+j}+\hat{a}_{l+j}^{\dagger}\hat{a}_{l}\right\rparen\cos\frac{n\pi}{L}l,\\\
&\sum_{l=1}^{L}\left\lparen\hat{a}_{l}^{\dagger}\hat{a}_{l+j}+\hat{a}_{l+j}^{\dagger}\hat{a}_{l}\right\rparen\sin\frac{n\pi}{L}l,\end{split}$
(16)

where $m\in\\{2,4,L\\}$, and $\left\lparen j,n\right\rparen\in\\{\left\lparen
1,2\right\rparen,\left\lparen 2,2\right\rparen,\left\lparen
1,L\right\rparen,\left\lparen 2,L\right\rparen\\}$. Note that these operators
are chosen so that they are not diagonal in the position or momentum basis.

Regarding the local reward $r_{t}$, we employ normalized LIOMs
${\\{\hat{I}_{n}^{+}/L\\}}_{1\leq n\leq n_{\mathrm{local}}}$, where we fix
$n_{\rm local}=4$ in the following. All $\hat{I}_{n}^{-}$ terms are excluded
since they are constantly zero not only for the initial and target states, but
also for any intermediate states $\rho_{t}$ evolved with the above unitaries.

###  V.2 Numerical results

Figure 4:  The learning curves for local preparation of the GGE. These figures
show the total rewards and the deviations of the LIOM densities at the end of
the episodes in evaluating the training progress. The horizontal axis shows
the number of parameter updates. (a) The total rewards. (b) The deviations of
the LIOM densities. The corresponding model is the transverse field Ising
model where the parameter is set so that the model is integrable. We added the
LIOMs $\hat{I}_{n}^{+}$ with separation $n=1,2,3,4$ to the reward of the RL.
The system size $L=120$. The lines show the median of the corresponding values
in evaluation with random seeds.

We now present the numerical results obtained by running the deep RL algorithm
to prepare prethermal pure quantum states that capture the characteristics of
the GGE. As we show the learning curve in Fig. 4(a), the deep RL agent
successfully learns to improve the local preparation protocol. This can be
more quantitatively understood from Fig. 4(b), which evaluates the absolute
difference in the expectation values of LIOMs between the target GGE and the
prepared state.

Figure 5:  The expectation values of the LIOMs in the protocol found by RL for
the GGE (blue). Each figure shows the LIOM with different separations. The red
horizontal lines show the corresponding expectation values of the target GGE.
The dotted black lines represent the relevant Gibbs ensemble, which shares the
same expectation values of total particle number and energy. The insets show
the absolute deviations in the expectation values
$\langle\hat{I}_{n}^{+}\rangle/L$ between the states in the protocol and the
target GGE. We added the LIOMs $\hat{I}_{n}^{+}$ with separation $n=1,2,3,4$
to the reward of the RL. Only LIOMs with separation $n\leq 4$ converge to the
expectation values of the GGE.

Let us further investigate the dynamics of the LIOMs generated by the
preparation protocol discovered by the RL agent. As shown in Fig. 5, the
behavior of LIOMs seem to be qualitatively different depending on $n$ in the
sense that, the only LIOMs with separation $n\leq n_{\mathrm{local}}=4$, which
we added to the reward, seem to converge to the expectation values of the GGE.

Figure 6:  The relaxation processes of the expectation values of correlation
functions (17), which are non-conserved quantities of $\hat{H}_{\mathrm{XX}}$
(12) (blue). Each figure shows the correlation function with different
separation $n$, which means the correlation function corresponds to the two
spins separated by $n-1$ sites. The red horizontal lines show the
corresponding expectation values of the target GGE, whereas the dotted black
lines represent the relevant Gibbs ensemble which shares the same expectation
values of total particle number and energy. Note that the system is controlled
with the learned protocol in $0\leq t\leq 40$, and subsequently undergoes free
evolution with $\hat{H}_{\mathrm{XX}}$. The only correlation functions with
separation $n\leq 4$ fluctuate around the expectation values of the GGE. Note
that the correlation function with $n=1$, where the Jordan-Wigner string does
not appear, is proportional to the LIOM with $n=1$, i.e., it is the conserved
quantity of $\hat{H}_{\mathrm{XX}}$.

Similar behavior can also be observed for non-conserved quantities such as the
correlation function

$\displaystyle\Gamma_{n}\coloneqq\frac{1}{L}\sum_{l}\left\lparen\hat{\sigma}_{l}^{+}\hat{\sigma}_{l+n}^{-}+\hat{\sigma}_{l+n}^{+}\hat{\sigma}_{l}^{-}\right\rparen,$
(17)

which can be seen as a representative of local operators acting on neighboring
$n+1$ sites. As we can see in Fig. 6, we see a notable convergence into the
GGE values for $n\leq n_{\rm local}$, while the errors seem to remain for
$n>n_{\rm local}$.

Figure 7:  Finite-size scaling of the distance function of the reduced density
operators between the prepared states and the target GGE. The dotted lines are
guides to eye that show power-law fitting $\bar{D}=\mathcal{O}(L^{-b})$, where
powers $b$ are summarized in Table 4 of Appendix. B. The error bars correspond
to the standard deviations of different protocols learned independently with
random seeds. The system is controlled by the deep RL agent at $0\leq t\leq
40$, and then undergoes free evolution by $\hat{H}_{\mathrm{XX}}$. Note that
the distance $\bar{D}$ concerns average over temporal, i.e., the free
evolution period, and the random training instance. The number of random seeds
is 5 at most.

Furthermore, we focus on the scaling of the error suppression. As we can see
from the finite-size scaling of the averaged distance function
$\bar{D}(\rho_{t},\rho_{\rm GGE})$ in Fig. 7, we observe that the distance of
the reduced density operators is suppressed as $\bar{D}=\mathcal{O}(L^{-b})$.
These behaviors are compatible with the scaling of the fluctuation of local
observables in the GME (6).

Here we conjecture that, while the errors for $LA\lesssim n_{\rm local}+1$
shall be suppressed polynomially even for larger system size, the errors for
$LA\gtrsim n_{\rm local}+1$ may saturate at finite values. This is because the
current local preparation scheme learns the LIOMs $\hat{I}_{n}^{\sigma}$ with
$n\leq n_{\rm local}$ to encode the prepared state into the “LIOM shell” only
for such conserved quantities. This means that, the prepared state fully
encodes the macroscopic properties of the tGGE, but not those of the GGE. We
numerically find that, the distance of the tGGE and GGE remains to be finite
even if the total system is in the thermodynamic limit [See Appendix F], which
is in agreement with Ref. [68] which investigated integrable region of the
transverse-field Ising chain. This supports our conjecture that the distance
between the GGE and the local-prepared state shall not be suppressed in the
thermodynamic limit. Meanwhile, it is possible that the error from the tGGE
itself is suppressed polynomially.

##  VI Conclusion

In this work, we have proposed a deep-RL-based quantum state preparation
framework for thermodynamic ensembles that solely relies on a few local
observables but not on global features such as the fidelity. The core idea is
to leverage the typicality of pure states in quantum many-body systems; the
macroscopic properties can be encoded simply via learning a few local
observables and undergoing free evolution. We have provided numerical
demonstrations that have successfully trained the deep RL agent to learn the
macroscopic properties of the Gibbs ensembles (Fig. 2) and the GGEs (Fig. 4).
We find that the accuracy of the prepared state has improved exponentially
with the system size for the former (Fig. 3) and polynomially for the latter
(Fig. 7), which is consistent with the argument of the typicality within a
given shell of local conserved quantities.

We envision three future directions of our work. The first important question
is the generalization of the local preparation protocol to include, e.g.,
dissipative terms, measurement and feedback, or postselection. We naturally
expect that the powerful explorability of the deep RL framework is not limited
to coherent control but could be applied to broader operation sets. Second, we
may consider the application of the local preparation protocol for the task of
Hamiltonian learning [71, 72, 73] by attempting to encode the macroscopic
properties using unitaries that do not explicitly contain the information of
Hamiltonian itself. Third, it is intriguing to seek how the local preparation
protocol is affected by various noises, such as the statistical noise
accompanied by sampling over observables. Efficient estimation methods such as
the randomized measurement schemes [74] shall be essential to boost the
training accuracy of the RL agent.

##  VII Acknowledgments

The authors wish to thank fruitful discussion with Ryusuke Hamazaki and Jiahao
Yao. S. B. is supported by Materials education program for the future leaders
in research, industry, and Technology (MERIT) of The University of Tokyo. N.Y.
wishes to thank JST PRESTO No. JPMJPR2119. Y.A. acknowledges support from the
Japan Society for the Promotion of Science through Grant Nos. JP19K23424 and
JP21K13859. T.S. is supported by JSPS KAKENHI Grant Number JP19H05796 and JST,
CREST Grant Number JPMJCR20C1, Japan. N.Y. and T.S. are also supported by
Institute of AI and Beyond of The University of Tokyo. We implemented the
exact time evolution in Sec. IV with QuSpin [75, 76] and the RL algorithm with
PyTorch [77]. The RL in Sec. V is performed on AI Bridging Cloud
Infrastructure (ABCI) of National Institute of Advanced Industrial Science and
Technology (AIST).

##  Appendix A Layout of the neural network and hyperparameters

Figure A1:  The abstract deep NN architecture used for the RL in the present work. At each time step $t$, the deep NN takes as input the observation $o_{t}$, the previous action $a_{t-1}$, and previous reward $r_{t-1}.$ After intermediate computations by fully-connected layers, LSTM, and the dueling network, the deep NN outputs the estimate of the Q function.  | Sec. IV | Sec. V
---|---|---
Linear + ReLU 1 | $240\rightarrow 512$ | $200\rightarrow 1024$
Linear + ReLU 2 | $512\rightarrow 512$ | $1024\rightarrow 1024$
LSTM | $(512,6,1)\rightarrow 512$ | $(1024,15,1)\rightarrow 1024$
Dueling network | $512\rightarrow 6$ | $1024\rightarrow 15$
Table 1:  The Input and Output sizes of each layer are shown in the former and latter of the arrow, respectively. The input of the NN consists of the action history, the one-hot representation of the previous action, and the previous reward. The input of Linear + ReLU 1 layer is only the action history. The one-hot representation of the previous action and previous reward is concatenated with the output of Linear + ReLU 2 before the LSTM layer.  Reward discount $\gamma$ |  | 0.997
---|---|---
Minibatch size |  | 324(Sec. IV)
|  | 380(Sec. V)
Sequence length |  | 40
Optimizer |  | Adam [78]
Optimizer setting | learning rate | $10^{-4}$
| $\varepsilon$ | $10^{-3}$
| $\beta$ | $(0.9,0.999)$
Replay ratio |  | 1
Gradient norms clip |  | 80
Table 2:  The hyperparameters used in the NNs. The agent performs updates on
batches of $(\text{minibatch size}\times\text{sequence length})$ observations.
Replay ratio means the effective number of times each experienced observation
is being replayed for the training. See Ref. [63] and its previous non-LSTM
version, Ref. [79], for the details of the hyperparameters. The other
parameters follow the ones in Ref. [62].

Here, we describe the architecture of the deep NN which is used to estimate
the action-value function (Q function). Figure A1 shows the overall picture;
the input for the LSTM at time step $t$ is the observation $o_{t}$, the
previous action $a_{t-1}$, and the previous reward $r_{t-1}$, whereas the
output is Q function whose optimal expression is given in Eq. (1) in the main
text. Refer to Table 1 for the size of the input and output of each layer and
Table 2 for the hyperparameter used in the NN. In the following, we further
describe the details of the structure.

The observation $o_{t}$ is chosen to be the action history $o_{t}=\left\lparen
a_{0},a_{1},\ldots,a_{t-1},-1,\ldots,-1\right\rparen$, which is fed to the
fully-connected layers. Fully-connected layers first perform linear
transformation, and then apply a non-linear activation function which is
chosen to be the rectified linear unit (ReLU) in the present work. The
intermediate output from the second fully-connected layer is concatenated with
the previous choice of action $a_{t-1}$ and the reward in the previous time
step $r_{t-1}$, and then fed to the Long Short-Term Memory (LSTM) layer.

The LSTM layer is introduced so that the network can refer to the history of
computational results at $t^{\prime}<t$ to estimate the Q function at time
step $t$. Namely, the input of the LSTM layer is not only the ones mentioned
above, but also its ”memory” including a _hidden state_ (short-term memory)
and a _cell state_ (long-term memory) [65]. Refer to literature such as Ref.
[80] for detailed information. This output memory is fed to the LSTM layer in
the next time step $t$, which enables the deep NN to successfully deal with
time-series inputs.

In the subsequent dueling network [81], the input is separated into two
branches. One branch evaluates the value of the observation $V(o)$, and the
other branch evaluates the advantage of actions regarding the observation
$A(o,a)=Q(o,a)-V(o)$. The output Q function of the dueling network is obtained
by summing the output of the two branches: $Q(o,a)=V(o)+A(o,a)$. This
separation may contribute to better training stability, faster convergence,
and better performance.

As a computational resource, we have used a single CPU and four GPUs (Intel
Xeon E5-2698 v4, $4\times$NVIDIA Tesla V100) in Sec. IV and two CPUs and four
GPUs ($2\times$Intel Xeon Gold 6148, $4\times$NVIDIA Tesla V100) in Sec. V,
respectively.

##  Appendix B Fitting parameters for the finite-size scaling of the
preparation accuracy

In Tables 3 and 4, we summarize the powers $b$ obtained by the fit for the
finite-size scaling of the local preparation accuracy in Sec. IV.2 and Sec.
V.2, respectively.

Subsystem size $LA$ | Power $b$
---|---
1 | 0.53(4)
2 | 0.54(4)
3 | 0.51(4)
4 | 0.48(4)
5 | 0.43(4)
Table 3:  The powers obtained by the fit for the result of the finite-size scaling of the distance function $\bar{D}$, given by Eq. (10), between the prepared states and the target Gibbs ensembles for the system size $L=10,12,14,16$ in Sec. IV.2.  subsystem size $LA$ | power $b$
---|---
1 | 0.50(3)
2 | 0.50(5)
3 | 0.46(4)
4 | 0.43(4)
5 | 0.41(4)
6 | 0.41(3)
7 | 0.41(3)
8 | 0.41(3)
9 | 0.41(3)
10 | 0.41(3)
Table 4:  The powers obtained by the fit for the result of the finite-size
scaling with the system size $L=60,72,84,96,108,120$ for the distance function
$\bar{D}$, given by Eq. (10), between the prepared states and the target GGE
in Sec. V.2.

##  Appendix C Control protocol discovered by the RL agent

Here, we briefly discuss the control protocol found by the RL agent. Figure C2
shows the control protocol whose learning curve shown in Fig. 2 in the main
text to realize a thermal pure state described by the Gibbs ensemble, which is
discussed in Sec. IV.2. It can be seen that the optimal protocol consists of
two stages: firstly choosing nontrivial actions that change the macroscopic
observables and then relaxing the system under free evolution using
$\hat{H}_{\rm Ising}$. While there is a small number of actions other than the
free relaxation at the late stage of the control protocol, we argue that they
do not contribute significantly to the result even if we substitute them with
$\hat{H}_{\rm Ising}$.

Figure C2:  The optimal action sequence $\\{a_{t}^{*}\\}$ for local
preparation of Gibbs ensemble discovered by the RL agent whose training
profile is provided in Fig. 2 in the main text. The six-fold action indices
correspond to 0: $\hat{H}_{\mathrm{Ising}}$, 1:
$\sum_{l}J\hat{\sigma}_{l}^{z}\hat{\sigma}_{l+1}^{z}+h\hat{\sigma}_{l}^{z}$,
2: $g\sum_{l}\hat{\sigma}_{l}^{x}$, 3: $\sum_{l}\hat{\sigma}_{l}^{y}$, 4:
$\sum_{l}\hat{\sigma}_{l}^{x}\hat{\sigma}_{l+1}^{y}+\hat{\sigma}_{l}^{y}\hat{\sigma}_{l+1}^{x}$,
5:
$\sum_{l}\hat{\sigma}_{l}^{y}\hat{\sigma}_{l+1}^{z}+\hat{\sigma}_{l}^{z}\hat{\sigma}_{l+1}^{y}$,
respectively. After quickly manipulating the macroscopic observables to encode
the property of the thermal ensemble, the agent mainly chooses to evolve the
system under the target Hamiltonian $\hat{H}_{\rm Ising}$.

##  Appendix D Relationship between energy shell width $\varepsilon_{\rm e}$
and scaling of preparation accuracy

Here we discuss the relationship between the energy shell width
$\varepsilon_{\mathrm{e}}$ and the scaling behavior of the accuracy of
prepared state. Recall that the distance function in the Eq. (10) of the main
text is given as

$\displaystyle D(\rho,\rho^{\prime})=\frac{\|\rho-\rho^{\prime}\|_{\rm
F}}{\sqrt{\|\rho\|_{\rm F}^{2}-\|\rho^{\prime}\|_{\rm F}^{2}}},$ (18)

where $\|A\|_{\rm F}=\sqrt{{\rm Tr}[A^{\dagger}A]}$ denoting the Frobenius
norm. In Fig. D3, we show how the scaling exponent $b$ defined from
$\bar{D}(\rho_{\rm Gibbs},\rho_{t})=O(d^{-b})$ varies according to
$\varepsilon_{\mathrm{e}}$. At $\varepsilon_{\mathrm{e}}=0.1$ for $L=10$, the
corresponding energy shell only includes a single eigenstate. Meanwhile,
$\varepsilon_{\mathrm{e}}=0.75$ corresponds to the extreme case where almost
all eigenstates below $\braket{\hat{H}_{\rm Ising}}_{\beta}$ are included in
the energy shell $[\braket{\hat{H}_{\rm
Ising}}_{\beta}-\varepsilon_{\mathrm{e}}L,\braket{\hat{H}_{\rm
Ising}}_{\beta}]$. In Sec. IV.2, we have chosen intermediate
$\varepsilon_{\mathrm{e}}$ so that the power $b$ is stable against the choice
of $\varepsilon_{\mathrm{e}}$.

Figure D3:  Relationship between the energy shell width $\varepsilon_{\rm e}$
and the scaling exponent $b$ of the distance $\bar{D}$, concerning the local
preparation of the Gibbs ensemble. The parameters of the system are identical
to those used in Sec. IV.2.

##  Appendix E Choosing Lagrange multipliers for the target GGE

The Lagrange multipliers for the GGE in Sec. V is chosen so that the
expectation values of local conserved quantities partly reproduce those of the
Gibbs ensemble, while some deviate from it. Concretely, we impose the
following equality:

$\displaystyle\operatorname{Tr}\left[\hat{I}_{n}^{+}\rho_{\mathrm{GGE,target}}\right]=\braket{\hat{I}_{n}^{+}}_{\beta}\
\ \ \,\text{for $n=0,~{}1$},$ (19)
$\displaystyle\operatorname{Tr}\left[\hat{I}_{n}^{+}\rho_{\mathrm{GGE,target}}\right]=\braket{\hat{I}_{n}^{+}}_{\beta}+\varepsilon_{\mathrm{I}}L\
\text{for $2\leq n<\frac{L}{2}$},$ (20)
$\displaystyle\operatorname{Tr}\left[\hat{I}_{n}^{-}\rho_{\mathrm{GGE,target}}\right]=\braket{\hat{I}_{n}^{-}}_{\beta},$
(21)

where ${\braket{\hat{I}_{n}^{\sigma}}_{\beta}}$ is the expectation values of
LIOMs of the Gibbs ensemble at the inverse temperature $\beta=0.4$. Note that
$\varepsilon_{\mathrm{I}}$ determines the deviation between the target GGE and
the Gibbs ensemble. For simplicity, we constantly take $\varepsilon_{\rm
I}=0.05$ for every $L$. In addition, we set
$\operatorname{Tr}\left[\hat{I}_{L/2}^{+}\rho_{\mathrm{GGE,target}}\right]=0$
because the fermionic tight-binding Hamiltonian obtained from the Jordan-
Wigner transformation (Eq. (13) in the main text) is anti-periodic when the
fermionic particle number $N_{\mathrm{f}}$ is even.

We remark that the LIOMs with $n=0,1$ correspond to the total particle number
and energy, respectively. Thus, this target GGE shares only the expectation
values of the total particle number and energy with the Gibbs ensemble.
Therefore, in order to prepare the subsystem whose size is larger than two, we
need to control additional LIOMs other than the particle number and energy.

##  Appendix F The distance between the GGE and the tGGE

Figure F4:  The distance, given by Eq. (10), between the target GGE in Sec.
V.2 and the corresponding tGGE with $n_{\mathrm{local}}=4$.

To quantify the difference between the tGGE and the target GGE considered in
Sec. V.2, we show the distance $D(\rho_{\rm GGE},\rho_{\rm tGGE})$ in Fig. F4
(see Eq. (10) for the definition). We observe that the tGGE and the GGE agree
well when $LA\leq n_{\mathrm{local}}+1$, whereas they deviate when
$LA>n_{\mathrm{local}}+1$. This result is compatible with Ref. [68], which
consider the integrable parameter region of transverse-field Ising chain.

## References

  * Nielsen and Chuang [2010] M. A. Nielsen and I. L. Chuang,  _Quantum computation and quantum information_ , 10th ed. (Cambridge University Press, Cambridge ; New York, 2010).
  * Ladd  _et al._ [2010] T. D. Ladd, F. Jelezko, R. Laflamme, Y. Nakamura, C. Monroe, and J. L. O’Brien, Quantum computers, Nature 464, 45 (2010).
  * Giovannetti  _et al._ [2011] V. Giovannetti, S. Lloyd, and L. Maccone, Advances in quantum metrology, Nature Photonics 5, 222 (2011).
  * Gisin  _et al._ [2002] N. Gisin, G. Ribordy, W. Tittel, and H. Zbinden, Quantum cryptography, Reviews of Modern Physics 74, 145 (2002).
  * Terhal and DiVincenzo [2000] B. M. Terhal and D. P. DiVincenzo, The problem of equilibration and the computation of correlation functions on a quantum computer, Physical Review A 61, 022301 (2000).
  * Poulin and Wocjan [2009] D. Poulin and P. Wocjan, Sampling from the Thermal Quantum Gibbs State and Evaluating Partition Functions with a Quantum Computer, Physical Review Letters 103, 220502 (2009).
  * Sugiura and Shimizu [2013] S. Sugiura and A. Shimizu, Canonical Thermal Pure Quantum State, Physical Review Letters 111, 010401 (2013).
  * [8] F. G. S. L. Brandao and M. J. Kastoryano, Finite correlation length implies efficient preparation of quantum thermal states, arXiv:1609.07877 .
  * Wu and Hsieh [2019] J. Wu and T. H. Hsieh, Variational Thermal Quantum Simulation via Thermofield Double States, Physical Review Letters 123, 220502 (2019).
  * [10] A. N. Chowdhury, G. H. Low, and N. Wiebe, A Variational Quantum Algorithm for Preparing Quantum Gibbs States, arXiv:2002.00055 .
  * Wang  _et al._ [2021] Y. Wang, G. Li, and X. Wang, Variational Quantum Gibbs State Preparation with a Truncated Taylor Series, Physical Review Applied 16, 054035 (2021).
  * Doria  _et al._ [2011] P. Doria, T. Calarco, and S. Montangero, Optimal Control Technique for Many-Body Quantum Dynamics, Physical Review Letters 106, 190501 (2011).
  * Caneva  _et al._ [2011] T. Caneva, T. Calarco, and S. Montangero, Chopped random-basis quantum optimization, Physical Review A 84, 022326 (2011).
  * Khaneja  _et al._ [2005] N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbrüggen, and S. J. Glaser, Optimal control of coupled spin dynamics: design of NMR pulse sequences by gradient ascent algorithms, Journal of Magnetic Resonance 172, 296 (2005).
  * Krotov [1996] V. F. Krotov,  _Global methods in optimal control theory_ , Monographs and textbooks in pure and applied mathematics No. 195 (M. Dekker, New York, 1996).
  * Mehta  _et al._ [2019] P. Mehta, M. Bukov, C.-H. Wang, A. G. Day, C. Richardson, C. K. Fisher, and D. J. Schwab, A high-bias, low-variance introduction to Machine Learning for physicists, Physics Reports 810, 1 (2019).
  * Carleo  _et al._ [2019] G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová, Machine learning and the physical sciences, Reviews of Modern Physics 91, 045002 (2019).
  * Yoshioka  _et al._ [2019] N. Yoshioka, Y. Akagi, and H. Katsura, Transforming generalized Ising models into Boltzmann machines, Physical Review E 99, 032113 (2019).
  * Nomura  _et al._ [2021] Y. Nomura, N. Yoshioka, and F. Nori, Purifying Deep Boltzmann Machines for Thermal Quantum States, Physical Review Letters 127, 060601 (2021).
  * Glasser  _et al._ [2018] I. Glasser, N. Pancotti, M. August, I. D. Rodriguez, and J. I. Cirac, Neural-Network Quantum States, String-Bond States, and Chiral Topological States, Physical Review X 8, 011006 (2018).
  * Gao and Duan [2017] X. Gao and L.-M. Duan, Efficient representation of quantum many-body states with deep neural networks, Nature Communications 8, 662 (2017).
  * Lu  _et al._ [2019] S. Lu, X. Gao, and L.-M. Duan, Efficient representation of topologically ordered states with restricted Boltzmann machines, Physical Review B 99, 155136 (2019).
  * Hibat-Allah  _et al._ [2020] M. Hibat-Allah, M. Ganahl, L. E. Hayward, R. G. Melko, and J. Carrasquilla, Recurrent neural network wave functions, Physical Review Research 2, 023358 (2020).
  * Sharir  _et al._ [2020] O. Sharir, Y. Levine, N. Wies, G. Carleo, and A. Shashua, Deep Autoregressive Models for the Efficient Variational Simulation of Many-Body Quantum Systems, Physical Review Letters 124, 020503 (2020).
  * Pfau  _et al._ [2020] D. Pfau, J. S. Spencer, A. G. d. G. Matthews, and W. M. C. Foulkes, Ab-Initio Solution of the Many-Electron schrödinger Equation with Deep Neural Networks, Physical Review Research 2, 033429 (2020).
  * Carleo and Troyer [2017] G. Carleo and M. Troyer, Solving the quantum many-body problem with artificial neural networks, Science 355, 602 (2017).
  * Lagaris  _et al._ [1997] I. E. Lagaris, A. Likas, and D. I. Fotiadis, Artificial neural network methods in quantum mechanics, Computer Physics Communications 104, 1 (1997).
  * Yoshioka and Hamazaki [2019] N. Yoshioka and R. Hamazaki, Constructing neural stationary states for open quantum many-body systems, Physical Review B 99, 214306 (2019).
  * Vicentini  _et al._ [2019] F. Vicentini, A. Biella, N. Regnault, and C. Ciuti, Variational Neural-Network Ansatz for Steady States in Open Quantum Systems, Physical Review Letters 122, 250503 (2019).
  * Hartmann and Carleo [2019] M. J. Hartmann and G. Carleo, Neural-Network Approach to Dissipative Quantum Many-Body Dynamics, Physical Review Letters 122, 250502 (2019).
  * Carrasquilla  _et al._ [2019] J. Carrasquilla, G. Torlai, R. G. Melko, and L. Aolita, Reconstructing quantum states with generative models, Nature Machine Intelligence 1, 155 (2019).
  * Torlai  _et al._ [2018] G. Torlai, G. Mazzola, J. Carrasquilla, M. Troyer, R. Melko, and G. Carleo, Neural-network quantum state tomography, Nature Physics 14, 447 (2018).
  * Torlai and Melko [2018] G. Torlai and R. G. Melko, Latent Space Purification via Neural Density Operators, Physical Review Letters 120, 240503 (2018).
  * Carrasquilla and Melko [2017] J. Carrasquilla and R. G. Melko, Machine learning phases of matter, Nature Physics 13, 431 (2017).
  * Yoshioka  _et al._ [2018] N. Yoshioka, Y. Akagi, and H. Katsura, Learning disordered topological phases by statistical recovery of symmetry, Phys. Rev. B 97, 205110 (2018).
  * Erdman and Noé [2022] P. A. Erdman and F. Noé, Identifying optimal cycles in quantum thermal machines with reinforcement-learning, npj Quantum Information 8, 1 (2022).
  * [37] P. A. Erdman and F. Noé, Driving black-box quantum thermal machines with optimal power/efficiency trade-offs using reinforcement learning, arXiv:2204.04785 .
  * Yoshioka  _et al._ [2021] N. Yoshioka, W. Mizukami, and F. Nori, Solving quasiparticle band spectra of real solids using neural-network quantum states, Communications Physics 4, 106 (2021).
  * Nagy and Savona [2019] A. Nagy and V. Savona, Variational quantum monte carlo method with a neural-network ansatz for open quantum systems, Phys. Rev. Lett. 122, 250501 (2019).
  * [40] S. Tibaldi, G. Magnifico, D. Vodola, and E. Ercolessi, Unsupervised and supervised learning of interacting topological phases from single-particle correlation functions, arXiv:2202.09281 .
  * Sutton and Barto [2018] R. S. Sutton and A. G. Barto,  _Reinforcement learning: an introduction_ , 2nd ed., Adaptive computation and machine learning series (The MIT Press, Cambridge, Massachusetts, 2018).
  * Bukov  _et al._ [2018] M. Bukov, A. G. R. Day, D. Sels, P. Weinberg, A. Polkovnikov, and P. Mehta, Reinforcement Learning in Different Phases of Quantum Control, Physical Review X 8, 031086 (2018).
  * Bukov [2018] M. Bukov, Reinforcement learning for autonomous preparation of Floquet-engineered states: Inverting the quantum Kapitza oscillator, Physical Review B 98, 224305 (2018).
  * Yao  _et al._ [2021] J. Yao, L. Lin, and M. Bukov, Reinforcement Learning for Many-Body Ground-State Preparation Inspired by Counterdiabatic Driving, Physical Review X 11, 031070 (2021).
  * Sivak  _et al._ [2022] V. V. Sivak, A. Eickbusch, H. Liu, B. Royer, I. Tsioutsios, and M. H. Devoret, Model-Free Quantum Control with Reinforcement Learning, Physical Review X 12, 011059 (2022).
  * Wang  _et al._ [2020] Z. T. Wang, Y. Ashida, and M. Ueda, Deep Reinforcement Learning Control of Quantum Cartpoles, Physical Review Letters 125, 100401 (2020).
  * Borah  _et al._ [2021] S. Borah, B. Sarma, M. Kewming, G. J. Milburn, and J. Twamley, Measurement-Based Feedback Quantum Control with Deep Reinforcement Learning for a Double-Well Nonlinear Potential, Physical Review Letters 127, 190403 (2021).
  * Niu  _et al._ [2019] M. Y. Niu, S. Boixo, V. N. Smelyanskiy, and H. Neven, Universal quantum control through deep reinforcement learning, npj Quantum Information 5, 33 (2019).
  * Ding  _et al._ [2021] Y. Ding, Y. Ban, J. D. Martín-Guerrero, E. Solano, J. Casanova, and X. Chen, Breaking adiabatic quantum control with deep learning, Physical Review A 103, L040401 (2021).
  * An  _et al._ [2021] Z. An, H.-J. Song, Q.-K. He, and D. L. Zhou, Quantum optimal control of multilevel dissipative quantum systems with reinforcement learning, Physical Review A 103, 012404 (2021).
  * Sgroi  _et al._ [2021] P. Sgroi, G. M. Palma, and M. Paternostro, Reinforcement Learning Approach to Nonequilibrium Quantum Thermodynamics, Physical Review Letters 126, 020601 (2021).
  * Fösel  _et al._ [2018] T. Fösel, P. Tighineanu, T. Weiss, and F. Marquardt, Reinforcement Learning with Neural Networks for Quantum Feedback, Physical Review X 8, 031084 (2018).
  * Rigol  _et al._ [2008] M. Rigol, V. Dunjko, and M. Olshanii, Thermalization and its mechanism for generic isolated quantum systems, Nature 452, 854 (2008).
  * Rigol  _et al._ [2007] M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Relaxation in a Completely Integrable Many-Body Quantum System: An Ab Initio Study of the Dynamics of the Highly Excited States of 1D Lattice Hard-Core Bosons, Physical Review Letters 98, 050405 (2007).
  * Vidmar and Rigol [2016] L. Vidmar and M. Rigol, Generalized Gibbs ensemble in integrable lattice models, Journal of Statistical Mechanics: Theory and Experiment 2016, 064007 (2016).
  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics, Advances in Physics 65, 239 (2016).
  * Goldstein  _et al._ [2006] S. Goldstein, J. L. Lebowitz, R. Tumulka, and N. Zanghì, Canonical Typicality, Physical Review Letters 96, 050403 (2006).
  * Popescu  _et al._ [2006] S. Popescu, A. J. Short, and A. Winter, Entanglement and the foundations of statistical mechanics, Nature Physics 2, 754 (2006).
  * Sugiura and Shimizu [2012] S. Sugiura and A. Shimizu, Thermal Pure Quantum States at Finite Temperature, Physical Review Letters 108, 240401 (2012).
  * [60] Y. Li, Deep Reinforcement Learning, arXiv:1810.06339 .
  * [61] P. Henderson, R. Islam, P. Bachman, J. Pineau, D. Precup, and D. Meger, Deep Reinforcement Learning that Matters, arXiv:1709.06560 .
  * [62] A. Stooke and P. Abbeel, rlpyt: A Research Code Base for Deep Reinforcement Learning in PyTorch, arXiv:1909.01500 .
  * Kapturowski  _et al._ [2019] S. Kapturowski, G. Ostrovski, J. Quan, R. Munos, and W. Dabney, Recurrent Experience Replay in Distributed Reinforcement Learning, in  _International Conference on Learning Representations_ (2019).
  * Mnih  _et al._ [2015] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski, S. Petersen, C. Beattie, A. Sadik, I. Antonoglou, H. King, D. Kumaran, D. Wierstra, S. Legg, and D. Hassabis, Human-level control through deep reinforcement learning, Nature 518, 529 (2015).
  * Hochreiter and Schmidhuber [1997] S. Hochreiter and J. Schmidhuber, Long Short-Term Memory, Neural Computation 9, 1735 (1997).
  * Note [1] Gibbs ensembles, which is our target, differs from microcanonical ensembles. Nevertheless, a microcanonical ensemble is locally equivalent to a Gibbs (canonical) ensemble for translation invariant short-ranged Hamiltonians, even if the system is finite [82, 83]. In other words, the microcanonical and the canonical expectation values almost coincide when we look at the subsystem that is not too large.
  * Cassidy  _et al._ [2011] A. C. Cassidy, C. W. Clark, and M. Rigol, Generalized Thermalization in an Integrable Lattice System, Physical Review Letters 106, 140405 (2011).
  * Fagotti and Essler [2013] M. Fagotti and F. H. L. Essler, Reduced density matrix after a quantum quench, Physical Review B 87, 245107 (2013).
  * Kim  _et al._ [2014] H. Kim, T. N. Ikeda, and D. A. Huse, Testing whether all eigenstates obey the eigenstate thermalization hypothesis, Physical Review E 90, 052105 (2014).
  * Note [2] Note that the chosen distance function can be efficiently computed for Slater determinants, which is used to express free-fermionic states obtained from Jordan-Wigner transformation. See Ref. [68] for the detailed property of this distance.
  * Rudinger and Joynt [2015] K. Rudinger and R. Joynt, Compressed sensing for Hamiltonian reconstruction, Physical Review A 92, 052322 (2015).
  * Bairey  _et al._ [2019] E. Bairey, I. Arad, and N. H. Lindner, Learning a Local Hamiltonian from Local Measurements, Physical Review Letters 122, 020504 (2019).
  * Anshu  _et al._ [2021] A. Anshu, S. Arunachalam, T. Kuwahara, and M. Soleimanifar, Sample-efficient learning of interacting quantum systems, Nature Physics 17, 931 (2021).
  * Huang  _et al._ [2020] H.-Y. Huang, R. Kueng, and J. Preskill, Predicting many properties of a quantum system from very few measurements, Nature Physics 16, 1050 (2020).
  * Weinberg and Bukov [2017] P. Weinberg and M. Bukov, QuSpin: a Python package for dynamics and exact diagonalisation of quantum many body systems part I: spin chains, SciPost Physics 2, 003 (2017).
  * Weinberg and Bukov [2019] P. Weinberg and M. Bukov, QuSpin: a Python package for dynamics and exact diagonalisation of quantum many body systems. Part II: bosons, fermions and higher spins, SciPost Physics 7, 020 (2019).
  * Paszke  _et al._ [2019] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, Pytorch: An imperative style, high-performance deep learning library, in  _Advances in Neural Information Processing Systems 32_, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019) pp. 8024–8035.
  * Kingma and Ba [2015] D. P. Kingma and J. Ba, Adam: A Method for Stochastic Optimization, in  _International Conference on Learning Representations_ (2015).
  * Horgan  _et al._ [2018] D. Horgan, J. Quan, D. Budden, G. Barth-Maron, M. Hessel, H. v. Hasselt, and D. Silver, Distributed Prioritized Experience Replay, in  _International Conference on Learning Representations_ (2018).
  * Sherstinsky [2020] A. Sherstinsky, Fundamentals of recurrent neural network (rnn) and long short-term memory (lstm) network, Physica D: Nonlinear Phenomena 404, 132306 (2020).
  * Wang  _et al._ [2016] Z. Wang, T. Schaul, M. Hessel, H. Hasselt, M. Lanctot, and N. Freitas, Dueling Network Architectures for Deep Reinforcement Learning, in  _Proceedings of The 33rd International Conference on Machine Learning_, Proceedings of Machine Learning Research, Vol. 48 (PMLR, New York, New York, USA, 2016) pp. 1995–2003.
  * [82] F. G. S. L. Brandao and M. Cramer, Equivalence of Statistical Mechanical Ensembles for Non-Critical Quantum Systems, arXiv:1502.03263 .
  * Tasaki [2018] H. Tasaki, On the local equivalence between the canonical and the microcanonical distributions for quantum spin systems, Journal of Statistical Physics 172, 905 (2018).