# Machine-learning-inspired quantum control in many-body dynamics

Meng-Yun Mao  College of Physics, Nanjing University of Aeronautics and
Astronautics, Nanjing 211106, China  Key Laboratory of Aerospace Information
Materials and Physics (NUAA), MIIT, Nanjing 211106, China  Zheng Cheng
College of Physics, Nanjing University of Aeronautics and Astronautics,
Nanjing 211106, China  Key Laboratory of Aerospace Information Materials and
Physics (NUAA), MIIT, Nanjing 211106, China  Liangsheng Li  National Key
Laboratory of Scattering and Radiation, Beijing 100854, China  Ning Wu
<EMAIL_ADDRESS>Center for Quantum Technology Research, School of Physics,
Beijing Institute of Technology, Beijing 100081, China  Key Laboratory of
Advanced Optoelectronic Quantum Architecture and Measurements (MOE), School of
Physics, Beijing Institute of Technology, Beijing 100081, China  Wen-Long You
<EMAIL_ADDRESS>College of Physics, Nanjing University of Aeronautics and
Astronautics, Nanjing 211106, China  Key Laboratory of Aerospace Information
Materials and Physics (NUAA), MIIT, Nanjing 211106, China

###### Abstract

Achieving precise preparation of quantum many-body states is crucial for the
practical implementation of quantum computation and quantum simulation.
However, the inherent challenges posed by unavoidable excitations at critical
points during quench processes necessitate careful design of control fields.
In this work, we introduce a promising and versatile dynamic control neural
network tailored to optimize control fields. We address the problem of
suppressing defect density and enhancing cat-state fidelity during the passage
across the critical point in the quantum Ising model. Our method facilitates
seamless transitions between different objective functions by adjusting the
optimization strategy. In comparison to gradient-based power-law quench
methods, our approach demonstrates significant advantages for both small
system sizes and long-term evolutions. We provide a detailed analysis of the
specific forms of control fields and summarize common features for
experimental implementation. Furthermore, numerical simulations demonstrate
the robustness of our proposal against random noise and spin number
fluctuations. The optimized defect density and cat-state fidelity exhibit a
transition at a critical ratio of the quench duration to the system size,
coinciding with the quantum speed limit for quantum evolution.

##  I Introduction

Nowadays it is imperative to achieve precise preparation of quantum many-body
states due to the urgent demand for applications in various fields, including
quantum computation [1, 2] and quantum simulation [3, 4, 5]. In particular,
quantum metrology, a rapidly advancing field of quantum information science,
is currently experiencing a surge of theoretical developments [6, 7, 8] and
experimental breakthroughs [9, 10, 11, 12]. A pivotal focus of a quantum
metrology scheme involves the preparation of optimal nonclassical states [13,
14]. Through the utilization of inherent quantum properties, such as
entanglement [15, 16, 17], coherence [18] and discord [19, 20], quantum
metrology can achieve the sensitivity for estimating unknown parameters that
surpasses the constraints imposed by the standard quantum limit of classical
strategies 111If there are $N$ qubits in separable states, it is known that
the uncertainty (that is, the inverse of the sensitivity) scales as ${\cal
O}(N^{-1/2})$, which is called the standard quantum limit, while quantum
physics allows the ultimate scalings being ${\cal O}(N^{-1})$, i.e., the
Heisenberg scaling, in the absence of noise and ${\cal O}(N^{-3/4})$ in the
presence of realistic decoherence. In fact, the enhanced sensitivity may even
approach the limits set by the Heisenberg uncertainty principle [22]. To this
end, various forms of entangled states have been generated in engineered many-
body systems to increase the phase sensitivity [23, 24]. However, a
comprehensive understanding of the relationship between quantum features and
ultimate scaling sensitivity beyond the standard quantum limit remains
elusive.

A widely employed method for the theoretical [25, 26, 27, 28, 29] and
experimental [30, 31, 32, 33, 34] preparation of quantum many-body states
involves utilizing unitary evolution. This process transforms the initial
state, typically the ground state of some simple Hamiltonian $\hat{H}_{i}$,
into the ground state of the target Hamiltonian $\hat{H}_{f}$. According to
the adiabatic theorem, the unitary transformation can be implemented with
arbitrary accuracy by changing the Hamiltonian sufficiently slowly [35, 36].
Nevertheless, the strict requirement of the adiabatic limit dictates
infinitely long evolution time, inevitably resulting in system decoherence and
the disappearance of entanglement [37]. In particular, crossing a quantum
critical point (QCP) in finite time challenges the adiabatic condition due to
the closing of the energy gap, which ultimately results in the formation of
excitations [38]. Consequently, a trade-off must be considered between the
evolution time and the excitations generated by the system. In this context,
tremendous efforts have recently been dedicated to achieving an optimal
passage through the non-adiabatic evolution, particularly across the QCP,
aiming to minimize unwanted excitations [39, 40] or to maximize fidelity [41,
42]. Shortcut to adiabaticity (STA) provides a way of finding fast
trajectories that connect the initial and final states by manipulating the
system’s parameters in a non-adiabatic fashion while still obtaining results
akin to those of an adiabatic process [43, 44, 45, 46, 47]. Given the
flexibility in selecting intermediate trajectories, the time-dependent control
parameters of a system can be adjusted in various STA protocols [48]. One
common method of STA involves introducing counterdiabatic driving to the
reference Hamiltonian that effectively compensates for the non-adiabatic
behavior. However, implementing this technique in quantum many-body systems
remains a challenging task for experimental execution.

Alternatively, a prevalent approach is quantum optimal control [49, 50, 51,
52, 53], where time-dependent control parameters of a system are fine-tuned
using optimal control theory. Governed by time-energy uncertainty relations,
the characteristic time scales during non-adiabatic evolution are encapsulated
by the quantum speed limit [54, 55, 56], which delves into the minimum time
required for quantum states to achieve specific predetermined objectives. This
is particularly crucial in the field of quantum information, where rapid
dynamics are often advantageous [57, 58, 56]. While the quantum optimal
control finds widespread applications in various systems, it is of fundamental
interest to formulate the control theory in a general framework. Over the past
few years, machine-learning technology has become an integral part of the
optimization theory and has been proved to be applicable to optimizing the
parameters of variational states in a variety of interacting quantum many-body
systems [59, 60, 61, 62, 63, 64]. In this work, we propose a promising and
generalizable dynamic control neural network (DCNN) to optimize the control
field passing through two fixed points of the quantum critical systems that
holds potential for testing and comparing approximate methods, as well as
bearing implications for condensed matter physics. To demonstrate the power of
our method, we focus on suppressing the defect density and improving the cat-
state fidelity during the passage across the critical point in the quantum
Ising model. Compared to the gradient-based power-law quench [65], the method
offers significant advantages over small system size and long-term evolution.
We also numerically demonstrate that our proposal is robust against the random
noise and the spin number fluctuations. The optimized defect density and cat-
state fidelity are uncovered to coincide with the quantum speed limit for
quantum evolution.

The rest of the paper is organized as follows. In Sec. II we present the
gradients of specific observables with respect to the time-dependent control
parameter in the many-body systems that can be transformed into free fermions.
Additionally, we introduce the DCNN method in this section. Section III is
devoted to the application of the DCNN to the suppression of the defect
density in the quantum Ising model. A thorough analysis of achieving the
optimized control fields is also discussed. The effectiveness of the protocol,
particularly under varying conditions such as random noise and fluctuations in
the number of spins, is evaluated through comprehensive numerical simulations.
Section IV delves into the improvement of the cat-state fidelity in the
quantum Ising model. Finally, the conclusion and outlook are given in Sec. V.

##  II Dynamical protocol in integrable many-body models

###  II.1 Quantum Ising-like models

The quantum Ising model has emerged as a prototypical model of quantum many-
body systems due to its analytical solvability. It thus serves as a valuable
tool for exploring emerging quantum phenomena [66, 67, 68, 69, 70, 71, 72] and
as a benchmark for evaluating the effectiveness of state-of-the-art approaches
and algorithms [73, 74, 75]. Furthermore, it is an especially compelling
option as a versatile platform for quantum simulation [76, 77, 78, 79]. The
exact simulation of Ising model has been implemented on quantum computers [80,
81]. Actually, the steps to build quantum circuits at a scale that could
provide an advantage for simulating the Ising Hamiltonian follow the same
strategy as the analytical solution of the models. Therefore, the method can
be extended to other integrable models. To this end, we focus on a class of
integrable many-body spin systems that can be mapped into free fermion models.
A variety of systems holding quantum phase transitions can be mapped into such
kinds of models, including the one-dimensional (1D) quantum Ising model and XY
model, which is regarded as one of two canonical quantum critical systems
[82], as well as the 1D and two-dimensional Kitaev models [83, 84], etc. These
cover a variety of spin models that are of interest in both quantum
information science and condensed matter physics.

Consider a closed quantum system which depends on time explicitly through a
control field $g(t)$ that can be varied arbitrarily with certain constraints.
The control field $g(t)$ is assumed to enter $\hat{H}[g(t)]$ via a term
$g(t)\hat{Y}$ with $\hat{Y}$ some time-independent operator. We consider a
quench starting from some fixed initial value $g_{i}=g(-T)$ and ending up with
a fixed final value $g_{f}=g(T)$ within the time interval $t\in[-T,T]$. In
this work, we mainly focus on quantum systems that themselves belong to, or
can be mapped to, $d$-dimensional free-fermion models. In general, this family
of arbitrary-dimensional time-dependent free-fermion Hamiltonians can be
written as

$\hat{H}(t)=\sum_{\vec{k}}\psi_{\vec{k}}^{\dagger}[\vec{d}_{\vec{k}}(g(t))\cdot\vec{\sigma}_{\vec{k}}]\psi_{\vec{k}},$
(1)

where
$\vec{\sigma}_{\vec{k}}=(\sigma_{\vec{k}}^{x},\sigma_{\vec{k}}^{y},\sigma_{\vec{k}}^{z})$
are the Pauli matrices acting on the mode $\vec{k}$ and
$\psi_{\vec{k}}=(a_{\vec{k}},b_{\vec{k}})^{T}$ are certain fermionic
operators. The function
$\vec{d}_{\vec{k}}(g)=(d_{\vec{k}}^{x}(g),d_{\vec{k}}^{y}(g),d_{\vec{k}}^{z}(g))$
is determined by specific models and the corresponding norm
$\varepsilon_{\vec{k}}=|\vec{d}_{\vec{k}}(g)|=\sqrt{d_{\vec{k}}^{x}(g)^{2}+d_{\vec{k}}^{y}(g)^{2}+d_{\vec{k}}^{z}(g)^{2}}$
gives the single-particle dispersion. Equation (1) can be diagonalized
directly within each $\vec{k}$-subspace as
$\hat{H}(t)=\sum_{\vec{k}}\hat{H}_{\vec{k}}(t)$ with

$\hat{H}_{\vec{k}}(t)=\varepsilon_{\vec{k}}(\Psi_{\vec{k}}^{\dagger}\Psi_{\vec{k}}-1),$
(2)

where $\Psi_{\vec{k}}=(A_{\vec{k}},B_{\vec{k}})^{T}$ is associated with
$\psi_{\vec{k}}$ by certain $\vec{k}$ dependent unitary transformation. Note
that both the dispersion $\varepsilon_{\vec{k}}$ and the quasi-fermion
operators $\Psi_{\vec{k}}$ are generally time-dependent through the control
field $g(t)$. The sub-ground state $|G_{\vec{k}}[g(t)]\rangle$ of mode
$\vec{k}$ satisfies

$A_{\vec{k}}|G_{\vec{k}}[g(t)]\rangle=B_{\vec{k}}|G_{\vec{k}}[g(t)]\rangle=0.$
(3)

The sub-excited state with the same fermion number parity as
$G_{\vec{k}}[g(t)]$ is obviously

$|\bar{G}_{\vec{k}}[g(t)]\rangle=A_{\vec{k}}^{\dagger}(t)B_{\vec{k}}^{\dagger}(t)|G_{\vec{k}}[g(t)]\rangle.$
(4)

We assume that the ground state is always nondegenerate in the process of the
quench and the energy gap is closed only at the QCP in the thermodynamic
limit. Here, both the global instantaneous ground state $|G[g(t)]\rangle$ and
time evolution operator $U(t,-T)$ are separable,

$\displaystyle|G[g(t)]\rangle$ $\displaystyle=$
$\displaystyle\prod_{\vec{k}}|G_{\vec{k}}[g(t)]\rangle,$ (5) $\displaystyle
U(t,-T)$ $\displaystyle=$ $\displaystyle\prod_{\vec{k}}U_{\vec{k}}(t,-T),$ (6)

where $U_{\vec{k}}(t,-T)=\mathcal{T}\exp(i\int_{-T}^{t}d\tau
H_{\vec{k}}(\tau))$ with $\mathcal{T}$ the time-ordered operator. The time-
evolved state is given by $|\phi(t)\rangle=U(t,-T)|G(g_{i})\rangle$, which is
also separable, i.e.,
$|\phi(t)\rangle=\prod_{\vec{k}}U_{\vec{k}}(t,-T)|G_{\vec{k}}(g_{i})\rangle$.

###  II.2 Controlling observables and their gradients.

Our aim is to minimize or maximize the final expectation value
$O(T)\equiv\langle\phi(T)|\hat{O}|\phi(T)\rangle$ of a general observable
$\hat{O}$. As will be discussed below, a key quantity utilized in the DCNN is
the gradient of $O(T)$ with respect to the control field $g(t)$ [85] (see
appendix A),

$\displaystyle\frac{\delta O(T)}{\delta g(t)}=2\imaginary\langle
G(g_{i})|\hat{O}(T)\hat{Y}(t)|G(g_{i})\rangle,$ (7)

where $\hat{A}(t)=U^{\dagger}(t,-T)\hat{A}U(t,-T)$ is the Heisenberg-picture
operator for any operator $\hat{A}$.

We further assume that the control term $\hat{Y}$ can be expressed as a
summation over even operators of independent modes

$\hat{Y}=\sum_{\vec{k}}\hat{Y}_{\vec{k}}.$ (8)

Here we are interested in the following two types of observables, i.e.,
$\hat{O}_{1}=\sum_{\vec{k}}\hat{O}_{1,\vec{k}}$ and
$\hat{O}_{2}=\prod_{\vec{k}}\hat{O}_{2,\vec{k}}$. The gradient of $O_{1}(T)$
has been calculated as (see Appendix A for detail)

$\displaystyle\frac{\delta O_{1}(T)}{\delta
g(t)}=2\imaginary\sum_{\vec{k}}\langle\phi_{\vec{k}}(T)|\hat{O}_{1,\vec{k}}|\bar{\phi}_{\vec{k}}(T)\rangle\langle\bar{\phi}_{\vec{k}}(t)|\hat{Y}_{\vec{k}}|\phi_{\vec{k}}(t)\rangle,$

where
$|\bar{\phi}_{\vec{k}}(t)\rangle=U_{\vec{k}}(t,-T)|\bar{G}_{\vec{k}}(g_{i})\rangle$
is the evolved sub-excited state of mode $\vec{k}$. The gradient of $O_{2}(T)$
can be calculated in a similar way,

$\displaystyle\frac{\delta O_{2}(T)}{\delta g(t)}$ $\displaystyle=$
$\displaystyle
2\prod_{\vec{p}}\langle\phi_{\vec{p}}(T)|\hat{O}_{2,\vec{p}}|\phi_{\vec{p}}(T)\rangle\times$
(10)
$\displaystyle\sum_{\vec{k}}\frac{\imaginary\langle\phi_{\vec{k}}(T)|\hat{O}_{2,\vec{k}}|\bar{\phi}_{\vec{k}}(T)\rangle\langle\bar{\phi}_{\vec{k}}(t)|\hat{Y}_{\vec{k}}|\phi_{\vec{k}}(t)\rangle}{\langle\phi_{\vec{k}}(T)|\hat{O}_{2,\vec{k}}|\phi_{\vec{k}}(T)\rangle}.\quad\quad$

###  II.3  Dynamic control neural network

We then introduce the DCNN method to optimize the observables in many-body
dynamics. As an extension of the control neural network 222Zheng Cheng, Meng-
Jiao Lyu, Takayuki Myo, Hisashi Horiuchi, Hiroshi Toki, Zhong-Zhou Ren,
Masahiro Isaka, Meng-Yun Mao, Hiroki Takemoto, Niu Wan, Wen-Long You, and Qing
Zhao, to appear, the DCNN is a neural network with an adaptive dynamic
structure and the flexible optimization strategy. This methodology not only
allows us to design more diverse control fields but also facilitates the
further minimization of specific observables, a capability not fully realized
by other methods, e.g., the direct gradient algorithm based on the power-law
quench [65]. To be concrete, the DCNN possesses a similar architecture to the
conventional neural network, which is composed of input, hidden, and output
layers. Two adjacent layers are connected by the weights, biases, and
activation function. The distinguishing features of the DCNN lie in its
dynamic structure of the network and the learning optimization strategy. The
dynamic structure of the network can be effectively executed by increasing the
number of units in each hidden layer. We opt for two hidden layers and
$\tanh(x)$ as the activation function.

Figure 1: The workflow of the DCNN and the visualization of the optimization
landscape. (a) The structure of the DCNN. We choose two hidden layers with
$n_{1}$ units in the first layer and $n_{2}$ in the second. $t$ is the input
of the neural network and $f(t)$ is the output. Each neuron in the hidden and
output layers has a corresponding weight and bias. The blue dots are the
original units in the neural network, the yellow dots are the newly increased
units and the dashed lines represent the corresponding increased connection.
(b) The optimization of the combined objective function $Q$ for the observable
$O$ and the loss $\mathcal{L}$. An initial trial form of the control field is
used as a starting point (blue triangle). $Q$ moves downhill (darker red
triangles) until the convergence is reached. To be noticed, the downhill
movement signifies a simultaneous decline in both $O$ and $\mathcal{L}$. (c)
The initial and final forms of the control field $g(t)$ with two fixed points
(black crosses).

The workflow of the DCNN is illustrated in Fig. 1. The neural network has two
hidden layers with $n_{1}$ units in the first layer and $n_{2}$ in the second.
The parameters of the neural network can be randomly initialized. For
simplicity, the parameters can be initialized by training the neural network
based on the empirical form of the control field. We then add $\delta n_{1}$
neurons to the first hidden layers and add $\delta n_{2}$ to the second. The
corresponding parameters of the newly added neurons should be randomly
initialized based on the average distribution. The input of the DCNN is $t$
and the output is the following specific function

$\displaystyle f(t)=W_{3}^{\rm{T}}\tanh.[W_{2}\tanh.(W_{1}t+B_{1})+B_{2}],$
(11)

where $W_{i}\ (i=1,2,3)$ are weight matrices, $B_{i}\ (i=1,2)$ are bias
matrices, and $\tanh.(x)$ is the function acting on each elements of the
matrix. The output $f(t)$ undergoes a linear transformation to produce the
control fields $g(t)$. Interestingly, increasing the number of neurons within
the neural-network-based approach exhibits a conceptual parallel in model
space expansion, akin to the methods employed in traditional optimal control
that augment the number of optimization parameters [87, 88]. Nevertheless, the
DCNN adaptively and flexibly explores model space, unlike the traditional
optimal control method that enhances it by requiring prior knowledge and
specific basis functions. The DCNN employs a pareto-optimization-based
strategy, adeptly balancing multiple objectives without compromising any
particular target. Due to the flexible optimization strategy of the DCNN, we
can impose desired constraints on the form of the control field. In our work,
the control field $g(t)$ has fixed values for initial moment $t=-T$ and final
moment $t=T$, which can be evaluated using the loss function

$\displaystyle\mathcal{L}=\Big{(}g(-T)-r\Big{)}^{2}+\Big{(}g(T)-s\Big{)}^{2},$
(12)

where $r$ and $s$ are the fixed points of $g(t)$ for $t=-T$ and $T$,
respectively. The other constraint, involving the extremum of the final
expectation value of a specific observable $O(T)$.

The DCNN is capable of conducting the optimization for multiple objectives,
aligning with approaches in quantum optimal control [89]. To be specific, the
classical approach to solve a multi-objective optimization problem is to
assign a weight $q_{i}$ to each normalized objective function so that the
problem is converted to a single objective problem with a scalar objective
function $Q$ [90]. In our approach, the objective function is not necessarily
normalized. We set $Q$ to be $Q=q_{1}O+q_{2}\mathcal{L}$, where $q_{1}$ and
$q_{2}$ are two arbitrary real numbers. Notably, the normalization of the
individual objective functions is not imperative, given the adaptability
provided by the adjustable coefficients $q_{i}$ ($i$ = 1 , 2). In order to
obtain the minimum value of $Q$, we adopt the gradient descent algorithm to
update all parameters,

$\displaystyle W_{i}$ $\displaystyle\leftarrow$ $\displaystyle
W_{i}+l^{W}_{i}\bigg{(}q_{1}\frac{1}{N}\frac{\partial O(T)}{\partial
W_{i}}+q_{2}\frac{\partial\mathcal{L}}{\partial W_{i}}\bigg{)},i=1,2,3,$
$\displaystyle B_{i}$ $\displaystyle\leftarrow$ $\displaystyle
B_{i}+l^{B}_{i}\bigg{(}q_{1}\frac{1}{N}\frac{\partial O(T)}{\partial
B_{i}}+q_{2}\frac{\partial\mathcal{L}}{\partial B_{i}}\bigg{)},i=1,2,$ (13)

for the next variation. Here, the learning rates $l^{W}_{i}$ and $l^{B}_{i}$
are adjustable parameters depending on the impact of the corresponding
parameters on $Q$.

The gradients of the $O(T)$ and $\mathcal{L}$ with respect to all the
parameters are given by

$\displaystyle\frac{\partial O(T)}{\partial W_{i}}$ $\displaystyle=$
$\displaystyle\int_{-T}^{T}dt\frac{\delta O(T)}{\delta g(t)}\frac{\partial
g(t)}{\partial W_{i}},$ $\displaystyle\frac{\partial O(T)}{\partial B_{i}}$
$\displaystyle=$ $\displaystyle\int_{-T}^{T}dt\frac{\delta O(T)}{\delta
g(t)}\frac{\partial g(t)}{\partial B_{i}},$ (14)

and

$\displaystyle\frac{\partial\mathcal{L}}{\partial
W_{i}}=2(g(-T)-r)\frac{\partial g(t)}{\partial
W_{i}}\bigg{|}_{t=-T}+2(g(T)-r)\frac{\partial g(t)}{\partial
W_{i}}\bigg{|}_{t=T},$ $\displaystyle\frac{\partial\mathcal{L}}{\partial
B_{i}}=2(g(-T)-r)\frac{\partial g(t)}{\partial
B_{i}}\bigg{|}_{t=-T}+2(g(T)-r)\frac{\partial g(t)}{\partial
B_{i}}\bigg{|}_{t=T}.$

The gradients of $g(t)$ with respect to its parameters are calculated using
the chain rule and are then normalized to mitigate the risk of gradient
explosion.

We only update the corresponding parameters of the newly increased neurons
according to Eq. (13), while keep the other parameters unchanged. Next, if
both the observable $O$ and the loss $\mathcal{L}$, calculated using the
updated parameters, are comparatively smaller than those obtained from the
original parameters, we will accept the new neurons as part of the neural
network, namely, $n_{1}\leftarrow n_{1}+{\delta n_{1}}$ and $n_{2}\leftarrow
n_{2}+{\delta n_{2}}$. Otherwise, the newly introduced units are reinitialized
randomly. We can iterate through the aforementioned steps to obtain the
desired control field. Adding more neurons into neural network architectures
involves directly with the efficiency and scalability of the model, presenting
numerical cost as an interesting aspect that warrants further study [91].
While it is challenging to provide a comprehensive overview due to the
complicate interplay between model complexity and computational resources, the
current insights suggest that computational time and memory requirements tend
to scale nearly linearly with the number of parameters (See more details in
appendix B). It is also important to note that, in order to prevent the neural
network from getting trapped in potential local minima, we can establish a
criterion for the neural network to adjust its architecture after an extended
period. Furthermore, we may approve the update even if the observable $O$ or
the loss $\mathcal{L}$, calculated using the updated parameters, is slightly
higher than those computed with the original parameters.

##  III Optimal suppression of defect generation in the quantum Ising model

###  III.1 Quantum Ising chain and gradient of the observable

In this following, we will concretely apply the DCNN method to physical
scenarios, specifically focusing on enhancing the suppression of defect
generation and improving the preparation of cat-states during the passage
across the QCP of the 1D quantum Ising model. For the sake of numerical
simplicity, we employ periodic boundary conditions for the spin chain [40,
43]. While certain local observables, i.e., local longitudinal magnetization,
may differ between periodic and open boundary conditions [92], it is expected
that the global observables, such as defect density and cat-state fidelity,
will demonstrate analogous behavior irrespective of the boundary conditions.
The Hamiltonian of a periodic quantum Ising chain with
$\vec{\sigma}_{N+1}=\vec{\sigma}_{1}$, is given by

$\displaystyle\hat{H}_{\mathrm{QIM}}=-\sum_{j=1}^{N}[\sigma_{j}^{x}\sigma_{j+1}^{x}+g(t)\sigma_{j}^{z}],$
(15)

where $g(t)$ is a global transverse magnetic field that we aim to control, and
for simplicity, we assume $N$ to be even. This model exhibits a quantum phase
transition at $g_{c}=1$ between the ferromagnetic phase for $0<g<1$ and the
paramagnetic phase for $g>1$.

The Hamiltonian $\hat{H}_{\mathrm{QIM}}$ can be diagonalized through a
standard procedure, which involves applying the Jordan-Wigner transformation
followed by the Fourier transformation on the fermionic operators. In this
study, we are interested in dynamical protocols starting from the paramagnetic
phase with $g_{i}>1$, ensuring that only the subspace with an even fermion
number parity is pertinent [93], yielding

$\displaystyle\hat{H}^{(+)}_{\mathrm{QIM}}$ $\displaystyle=$
$\displaystyle\sum_{k\in K}\hat{H}_{k},$ $\displaystyle\hat{H}_{k}$
$\displaystyle=$ $\displaystyle
2(c^{\dagger}_{k},c_{-k})\left(\begin{array}[]{cc}g+\cos k&\sin k\\\ \sin
k&-g-\cos k\\\ \end{array}\right)\left(\begin{array}[]{c}c_{k}\\\
c^{\dagger}_{-k}\\\ \end{array}\right),$ (20)

where $K=\\{\pi/N,3\pi/N,\cdots,\pi-\pi/N\\}$ and the $c_{k}$’s are fermionic
operators satisfying $\\{c_{k},c_{k^{\prime}}\\}=0$ and
$\\{c_{k},c^{\dagger}_{k^{\prime}}\\}=\delta_{kk^{\prime}}$. Note that for
open boundary conditions, the necessary adjustment to the momentum set $K$
involves solving the corresponding transcendental equation [94]. The ground
and excitation states of the mode Hamiltonian $\hat{H}_{k}$ are given by

$\displaystyle|G_{k}\rangle$ $\displaystyle=$
$\displaystyle\cos\frac{\theta_{k}}{2}|\mathrm{vac}\rangle_{k}+\sin\frac{\theta_{k}}{2}|k,-k\rangle,$
$\displaystyle|\bar{G}_{k}\rangle$ $\displaystyle=$
$\displaystyle\sin\frac{\theta_{k}}{2}|\mathrm{vac}\rangle_{k}-\cos\frac{\theta_{k}}{2}|k,-k\rangle,$
(22)

where $|\mathrm{vac}\rangle_{k}$ is the vacuum state of $c_{\pm k}$ and
$|k,-k\rangle=c^{\dagger}_{k}c^{\dagger}_{-k}|\mathrm{vac}\rangle_{k}$ . The
angle $\theta_{k}$ is determined by $\tan\theta_{k}=-\sin k/(g+\cos k)$, where
$\Lambda_{k}(t)=2\sqrt{g^{2}(t)+2g(t)\cos k+1}$ is the single-particle
dispersion with which one has
$\hat{H}_{k}|G_{k}\rangle=-\Lambda_{k}|G_{k}\rangle$ and
$\hat{H}_{k}|\bar{G}_{k}\rangle=\Lambda_{k}|\bar{G}_{k}\rangle$. In the
present problem, the operator $\hat{Y}=-\sum_{j}\sigma^{z}_{j}$ gives

$\displaystyle\hat{Y}_{k}=-2\sum_{k\in
K}(c^{\dagger}_{k}c_{k}+c^{\dagger}_{-k}c_{-k}-1).$ (23)

Defects or domain walls will be inevitably produced during the transition from
the paramagnetic phase to the ferromagnetic phase. The number of defects is
measured by the operator

$\displaystyle\hat{D}=\frac{1}{2}\sum_{j=1}^{N}(1-\sigma_{j}^{x}\sigma_{j+1}^{x})=\sum_{k>0}\hat{D}_{k},$
(24)

where $\hat{D}_{k}$ is represented by the corresponding matrix form

$\mathcal{D}_{k}=1+\left(\begin{array}[]{cc}\cos k&-\sin k\\\ -\sin k&-\cos
k\end{array}\right)$ (25)

in the two-dimensional space spanned by
$\\{|\mathrm{vac}\rangle_{k},|k,-k\rangle\\}$. According to Eq. (II.2), the
gradient of $D(T)$ with respect to $g(t)$ is

$\displaystyle\frac{\delta D(T)}{\delta g(t)}$ $\displaystyle=$ $\displaystyle
4\textrm{Im}\sum_{k>0,k\in K}\langle\phi_{k}(T)|\cos k\sigma_{k}^{z}-\sin
k\sigma_{k}^{x}|\bar{\phi}_{k}(T)\rangle$ (26)
$\displaystyle\langle\bar{\phi}_{k}(t)|\sigma_{z}^{k}|\phi_{k}(t)\rangle.$

###  III.2 DCNN guided by the final defect density

We consider a quench from ${g(-T)}=2$ to ${g(T)}=0$ within the time interval
$[-T,T]$. We focus on a power-law quench characterized by the function
$g(t)=1-|t/T|^{r}\textrm{sgn}(t)$, which has been demonstrated to effectively
reduce the final defect density, defined as $\rho(T)\equiv D(T)/N$ [65].
Notably, this power-law quench strategy surpasses the conventional linear
quench [40], the assisted adiabatic passage [43] and the local adiabatic
evolution approach [95] in minimizing $\rho(T)$. To test the DCNN method, we
firstly consider the case of $T=50$ and $N=50$. In the initial neural network,
there are two hidden layers with $n_{1}=10$ neurons in the first layer and
$n_{2}=25$ neurons in the second layer. To improve the learning efficiency of
the neural network, we employ pre-training based on a gradient-oriented power-
law quench to establish the initial parameters. We then set ${\delta n_{1}}=5$
and ${\delta n_{2}}=10$. These values are adaptable and can be increased
throughout the learning process. The parameters for the newly added neurons in
these layers are randomly initialized, following a uniform distribution
ranging from $-0.1$ to $0.1$. We set the learning rates to be
$l^{W}_{i}=l^{B}_{i}=0.05$ for $\ i=1,2$ and $l^{W}_{3}=0.5$. Due to the
significant impact of $W_{3}$ on the $Q$, the learning rate of $W_{3}$ is set
at a value tenfold higher than that of the other parameters. The weights
assigned to $Q$ are $q_{1}=5$ and $q_{2}=1$, and we ensure the normalization
of gradients throughout the learning process.

The learning results of the DCNN for $T=0.5,N=50$ and $T=50,N=50$ are shown in
Fig. 2. The number of neurons in the first layer $n_{1}$ exhibits a stepwise
increase during the learning process. Each step signifies a transformation in
the neural network’s scale, leading to an enlargement of the model’s
representational space. During these steps, new nodes are introduced with
different initializations. The primary objective of this process is to
incrementally decrease the value of $\rho(T)$ through successive modifications
of $g(t)$ facilitated by these additional nodes. The interplay between
reducing the final defect density and the system’s evolution through two fixed
points presents an intriguing scenario. In certain cases, reducing $\rho(T)$
may result in the control field $g(t)$ diverging more substantially from these
fixed points. To accelerate the convergence rate of the neural networks, we
will selectively abandon the optimized form of the control field, although the
trade-off could potentially become less pronounced as the neural network’s
learning process advances. In this context, the DCNN exhibits a pronounced
ability to navigate multiple objectives and thereby effectively mitigates
oscillatory behaviors between competing objectives, fostering a more rapid
convergence. In Fig. 2, we observe an oscillatory decrease in $\rho(T)$ with
increasing learning iteration number $\mathcal{N}$, indicating significant
modifications to the neural network due to the introduction of additional
modes. Note that outliers resulting from significant modifications in the
learning process should not be a concern, as the optimization strategy
automatically discards these values. As the learning process progresses in
both scenarios, $\rho(T)$ eventually exhibits convergence with minimal
fluctuations. Specifically, $\rho(T)$ converges to a finite value for $T=0.5$,
while for $T=50$ it stabilizes around a minimal value.

###  III.3 Numerical results and discussions

In Fig. 3, we present the optimized control field $g(t)$ obtained by both the
DCNN method and the gradient-based power-law quench. In contrast to the power-
law control, the resulting $g(t)$ lacks symmetric structure and passes through
the QCP at a specific time $t^{*}(<0)$. Interestingly, we observe that $t^{*}$
tends to increase as $T$ increases. For the long-term evolution, it remains
imperative to navigate through the QCP as slowly as feasible. Nevertheless, in
both instances, the slowest rate of change in $g(t)$ occurs at $t=0$. These
characteristics of $g(t)$ could provide insights into the optimal timing and
approach for traversing the critical point, thereby minimizing defect
generation in practical pulse design applications. Furthermore, such insights
may hold the potential for generalization to state preparation protocols
across various quantum critical systems.

Figure 2: The learning results of the DCNN for (a) $T=0.5,N=50$ and (b)
$T=50,N=50$. The number of neurons in the first hidden layer $n_{1}$ (orange
dash-dot line) illustrates the change of the structure of the neural network
with the learning iterations $\mathcal{N}$. The blue line shows the specific
optimization process of the final defect density $\rho(T)$.  Figure 3: The
comparison between the optimized forms of the control field $g(t)$ based on
the DCNN (red solid line) and the gradient-based power-law quench (purple
dash-dot line) for (a) $T=0.5,N=50$ and (b) $T=50,N=50$. In both panels, the
horizontal grey dashed lines indicate QCPs, while the grey dashed vertical
lines identify $t^{*}$.

The optimized final defect density $\rho(T)$ using the DCNN is plotted in Fig.
4(a) for various system sizes with $T/N$ ranging from 0.01 to 1. To assess the
relative effectiveness of different optimization methods, we plot in Fig. 4(b)
the ratio $R_{\rho}$ of the defect density optimized by the DCNN to that
achieved through the gradient-based power-law quench. We see that the DCNN
outperforms the optimal power-law quench across a broad range of parameters,
particularly for small-sized systems and in the long-duration limit. Similar
to the case of the optimal power-law profile, we observe a crossover in the
scaled time $T/N$ beyond which the optimized defect density undergoes a sharp
decline. The simultaneous occurrence of this crossover for different values of
$T$ and $N$ suggests that the crossover time $T_{\rm c}$ is proportional to
$N$, contrasting with the Kibble-Zurek scale for the linear quench, where
$T_{\rm KZ}$ scales as $N^{2}$ [40]. This occurrence can be ascribed to the
implications of the quantum speed limit, which establishes a connection
between the maximum speed of evolution and the system’s energy uncertainty and
mean energy. Prior studies suggest that the quantum evolution in a many-body
model primarily involves two lowest states in the scenario of small
excitations [96, 97, 98]. This enables effective mapping into two-level
Hamiltonians [99, 100, 101, 102]. In the even subspace of $k$ mode in the 1D
transverse Ising model, the mode Hamiltonian $\hat{H}_{k}$ in Eq.
(LABEL:eq:Hk) can be reformulated in a Landau-Zener-type form as
$\hat{H}_{k}=\gamma_{k}\sigma_{k}^{z}+\omega_{k}\sigma_{k}^{x}$ with
$\gamma_{k}=2(g+\cos k)$ and $\omega_{k}=-2\cos k$. There exists an intrinsic
quantum speed limit $T_{\rm QSL}$ to drive a general initial state to a final
desired state [51], expressed by

$\displaystyle\tan 2\omega_{k}T_{\rm QSL}\propto\frac{1}{\omega_{k}}.$

For the quench with $t\in[-T,T]$ in the present scenario, the lowest mode
$k_{N}$ specified by $\omega_{k_{N}}\approx-2\pi/N$ yields
$2\omega_{k_{N}}T_{\rm QSL}\approx-\pi/2$ or $T_{\rm QSL}/N\approx 1/8$ for
sufficiently large $N$, serving as a lower bound for $T_{\rm QSL}$ for the
entire system, considering all modes $k$.

Figure 4: (a) The final defect density $\rho(T)$ with respect to different $T$
and $N$ under the optimization by the DCNN. (b) The ratio $R_{\rho}$ of the
final defect density optimized by the DCNN relative to that obtained by the
gradient-based power-law quench method. Figure 5: (a)The optimized control
fields $g(t)$ obtained by the DCNN protocol with different initial values
$g_{i}$ for $T=50,N=50$. Inset magnifies the control fields around $t^{*}$.
(b) The absolute value of the relative offset of $g(t)$ to $g_{c}$ (dashed
lines) and the energy gap $\Delta E$ (solid lines with markers) with respect
to $t$. Inset shows the corresponding final defect density $\rho(T)$ for
different $g_{i}$.

The gradient-based power-law quench method, which depends on a specific power-
law scheme to pass through fixed points, faces limitations in adaptability
across varied initial and final state scenarios. This limitation is
particularly pronounced when these states do not exhibit symmetry relative to
the critical point, thereby complicating the determination of a suitable
power-law function. In contrast, the DCNN offers a more adaptable framework,
enabling more straightforward modifications to its optimization strategy. This
adaptability suggests that the DCNN approach may provide more nuanced and
flexible solutions compared to the gradient-based power-law quench method,
especially regarding its application across diverse quenching scenarios. We
plot the optimized control fields with different initial values of $g_{i}$ in
Fig. 5 (a) for $T=50$ and $N=50$. It is evident that $g(t)$ smoothly passes
through the QCP at some $t^{*}$ around $t=0$ at a relatively slow pace. As
$g_{i}$ increases while keeping $g_{f}$ constant, mitigating non-adiabatic
effects within a finite duration becomes challenging. The optimized control
fields strikes a balance between a deliberate traversal of the QCP and
avoiding excessive speed in other regions [39]. To investigate the role of the
lowest mode in the optimization process, we plot the absolute value of the
relative offset of $g(t)$ to $g_{c}$ and the energy gap $\Delta E$ in Fig. 5
(b). One finds a discernible correlation between the deviation of $g(t)$ from
$g_{c}$ and the energy gap associated with the lowest mode, encapsulated by
the approximate relation $|g(t)-g_{c}|\approx 0.25\Delta E$.

To explore potential applications, we assess the robustness of our protocol to
spin number fluctuations [50], considering cases with $T=50,N=50$ and
$T=0.5,N=50$ as examples. In experimental implementations, it is inevitable to
encounter discrepancies in the actual number of spins present in the system.
We thus introduce a number fluctuation parameter, $\delta N=\pm 4$. The final
defect density remains below $9.13\times 10^{-5}$ for $T=50,N=50$ and below
$0.1798$ for $T=0.5,N=50$, indicating the robustness of our proposal against
spin number fluctuations. To delve deeper into the effects of random noise, we
examine the additive white Gaussian noise model, denoted as
$\mathcal{A}_{G}[g(t),\textrm{SNR}]$. We incorporate it into the control field
using the relation:

$g_{q}(t)=g(t)+\mathcal{A}_{G}[g(t),\textrm{SNR}],$ (27)

where $q$ represents each random generator
$\mathcal{A}_{G}[g(t),\textrm{SNR}]$, with
$\textrm{SNR}=10\log_{10}(P_{s}/P_{n})$, where $P_{s}$ and $P_{n}$ are the
power of signal and noise, respectively. We evaluate our protocol’s robustness
against white Gaussian noise by analyzing the mean relative offsets,
$\delta\rho$, from the ideal defect densities, conducting $\mathcal{R}$
numerical simulations with $\textrm{SNR}=10$ for $T=0.5,N=50$ and $T=50,N=50$.
Fig. 6 depicts $\delta\rho$ versus $\mathcal{R}=50p$, where $p$ is the
simulation index, showing convergence of defect densities to $3.57\times
10^{-5}$ and 0.1798 for the respective cases. These results highlight our
scheme exhibits a high degree of robustness against random noise.

Figure 6:  Mean relative offsets $\delta\rho$ of the final defect density from
the ideal value subjected to white Gaussian noise, as a function of the total
number of numerical simulations $\mathcal{R}$, with $\textrm{SNR}=10$. (a)
$T=0.5$, $N=50$, and (b) $T=50$, $N=50$. Here, we define $\mathcal{R}=50p$,
where $p=155$ for (a) and $p=438$ for (b). Dashed lines in both panels
indicate the convergent values.

##  IV Optimal cat-state preparation in the quantum Ising model

In this section, we employ the DCNN to tackle a more complex task: optimally
preparing a cat state in the quantum Ising chain. The Schrödinger cat states
have attracted much attention from researchers due to their fundamental
implications. The quantum superpositions of two macroscopically distinct
states are not only interesting for testing the fundamentals of quantum
mechanics such as the quantum-to-classical transition [103, 104], but also a
valuable resource for quantum information processing [105, 106, 107].
Noteworthy examples of maximally entangled cat states include the NOON state
[108] and Greenberger-Horne-Zeilinger (GHZ) state [109], which can push the
estimation precision to the Heisenberg limit [110, 111, 112]. Of particular
interest are spin cat states for use in atomic clocks [113] and magnetometers
[114] due to the potential applications in diverse fields from materials
science to medical diagnostics [115]. However, the realization of cat states
in practical experiments faces challenges due to their sensitivity to
decoherence and the delicate entanglement properties [116].

The target state we aim to achieve is

$\displaystyle|\mathrm{cat}\rangle=\frac{1}{\sqrt{2}}(|\rightarrow\cdots\rightarrow\rangle+|\leftarrow\cdots\leftarrow\rangle),$
(28)

which is an equally weighted superposition of the two degenerate fully
polarized states $|\rightarrow\cdots\rightarrow\rangle$ and
$|\leftarrow\cdots\leftarrow\rangle$. Our goal is to identify an optimized
function $g(t)$ that maximizes the final fidelity, defined as

$\displaystyle F(T)=|\langle\phi(T)|\mathrm{cat}\rangle|^{2}.$ (29)

It is important to note that maximizing this fidelity poses a more stringent
requirement than simply minimizing the final defect density. This is because
any linear combination of $|\rightarrow\cdots\rightarrow\rangle$ and
$|\leftarrow\cdots\leftarrow\rangle$ would result in vanishing defects.

It is shown in Ref. [93] that $|\mathrm{cat}\rangle$ can be written in the
momentum space as a product state

$\displaystyle|\mathrm{cat}\rangle=\prod_{k\in K}|\chi_{k}\rangle$ (30)

with
$|\chi_{k}\rangle\equiv\sin\frac{k}{2}|\mathrm{vac}\rangle_{k}+\cos\frac{k}{2}|k,-k\rangle$.
Therefore, $F(T)$ is actually the expectation value of the operator
$\hat{F}\equiv\prod_{k\in K}|\chi_{k}\rangle\langle\chi_{k}|$ in the final
state $|\phi(T)\rangle$. From Eq. (10), the dependence of $F(T)$ with respect
to $g(t)$ is delineated as

$\displaystyle\frac{\delta F(T)}{\delta g(t)}$ $\displaystyle=$ $\displaystyle
2\textrm{Im}\prod_{k}|\langle\phi_{k}(T)|\chi_{k}\rangle|^{2}$ (31)
$\displaystyle\sum_{k^{\prime}}\frac{\langle\chi_{k}|\bar{\phi}_{k^{\prime}}(T)\rangle\langle\bar{\phi}_{k^{\prime}}(t)|\hat{Y}_{k^{\prime}}|\phi_{k^{\prime}}(t)\rangle}{\langle\chi_{k}|\phi_{k^{\prime}}(T)\rangle}.$

Within the DCNN framework, we also examine the quench from $g(-T)=2$ to
$g(T)=0$ over the time interval $[-T,T]$. The objective function is
$Q=q_{1}F(T)+q_{2}\mathcal{L}$. The update of all parameters given by Eq. (13)
will undergo corresponding adjustments to reflect changes in $Q$. The
optimized fidelity results $F(T)$ obtained through the DCNN for sizes
$N=24,50,100$ and the ratio $T/N$ ranging from 0.01 to 1 are presented in Fig.
7. As expected, the achieved fidelity increases as the duration $T$ increases
and approaches unity for larger $T/N$. In the intermediate quench time range,
specifically between $0.14<T/N<0.18$, a notable surge in $F(T)$ is observed,
which is accompanied by a marked decrease in final defect density, as
illustrated in Fig. 4. It is noteworthy that this rise in $F(T)$ becomes more
pronounced with the increase in system size. For shorter durations, the
fidelity always remains extremely low, a consequence of the inherent
constraints imposed by the quantum speed limit.

Fig. 8 displays the optimized control fields $g(t)$, obtained either by
maximizing cat-state fidelity $F(T)$ or by minimizing final defect density
$\rho(T)$. In Fig. 8(a), a notable difference is observed between the two
control fields for short durations. Specifically, the $g(t)$ that minimizes
final defect density leads to considerably lower fidelity. Conversely, for
longer durations, the two profiles are nearly identical in Fig. 8(b), with
each leading to either very low final defect density or significantly high
fidelity. This observation implies a complex relationship between final defect
density and state fidelity, indicating that they are not straightforwardly
correlated.

Figure 7: The fidelity $F(T)$ of the cat state is evaluated under various
values of $T$ and $N$, optimized by the DCNN.  Figure 8: Comparison of the
control fields $g(t)$ optimized for the cat-state fidelity $F(T)$ (purple
dash-dot line) and the final defect density $\rho(T)$ (red solid lines) with
(a) $T=0.5,N=50$ and (b) $T=50,N=50$. In both panels, the horizontal grey
dashed lines indicate QCPs, while the grey dashed vertical lines identify
$t^{*}$ of $g(t)$ guided by the cat-state fidelity.

##  V Conclusion and outlook

In this study, we introduce a dynamic control neural network (DCNN) approach
to suppress the final defect density during a passage across a quantum
critical point (QCP). Utilizing the dynamic architecture of the neural network
combined with a strategic control mechanism, we successfully generate control
fields that navigate through two fixed points while optimizing the desired
objective. By considering various evolution time $T$ and system sizes $N$, our
method achieves a lower final defect density compared to traditional gradient-
based power-law quench methods. Meanwhile, our study reveals that the
optimized final defect density undergoes a crossover at a specific critical
ratio of quench duration to system size. This transition closely aligns with
the quantum speed limit for evolution, which is profoundly influenced by the
lowest mode. We numerically demonstrate that our protocol is robust against
the random noise and the spin number fluctuations. Additionally, we
demonstrate the effectiveness of our method using control fields that are
optimized across a range of initial values. By analyzing variational control
field forms and comparing them to power-law quench methods, we offer insights
for experimental pulse design. The versatility of the DCNN is further
highlighted by its application to enhancing cat-state fidelity, providing
similar optimized solutions indicative of the quantum speed limit. Our
findings reveal a non-one-to-one correspondence between final defect density
and fidelity. To attain high-quality quantum many-body states, it is crucial
to exceed the critical ratio of quench duration to system size.

Noting that the DCNN method introduced in this study extends beyond the
preparation of spin cat states, holding promise for a broad range of
applications in quantum technologies. With the rapid advancements in quantum
computation [117, 118, 119, 120] and quantum metrology [121, 122, 123, 124],
the DCNN’s capability for optimizing complex quantum systems positions it as a
valuable tool for enhancing precision in quantum information science. The
application of the DCNN method to multi-objective optimization in quantum
control demonstrates its potential applicability, especially in complex
scenarios such as quantum many-body systems and open quantum systems.

###### Acknowledgements.

This work is supported by the National Natural Science Foundation of China
(NSFC) under Grant No. 12174194 and stable supports for basic institute
research under Grant No. 190101. N. W. was supported by the National Key
Research and Development Program of China under Grant No. 2021YFA1400803.

##  Appendix A The gradients of $O(T)$ with respect to $g(t)$

The gradients of the final expectation value $O(T)$ with respect to the
control field $g(t)$ can be written as

$\displaystyle\frac{\delta O(T)}{\delta g(t)}=2\real\langle
G(g_{i})|U^{\dagger}(T,-T)\hat{O}\frac{\delta U(T,-T)}{\delta
g(t)}|G(g_{i})\rangle,\quad$ (32)

where

$\displaystyle\frac{\delta U(T,-T)}{\delta g(t)}$ $\displaystyle=$
$\displaystyle-iU(T,-T)U^{\dagger}(t,-T)\frac{\partial H(t)}{\partial
g(t)}U(t,-T)$ (33) $\displaystyle=$ $\displaystyle-
iU(T,-T)U^{\dagger}(t,-T)\hat{Y}U(t,-T).$

Combining Eq. (32) with Eq. (33), we derive the following expression:

$\displaystyle\frac{\delta O(T)}{\delta g(t)}=2\imaginary\langle
G(g_{i})|\hat{O}(T)\hat{Y}(t)|G(g_{i})\rangle.$ (34)

For simplicity, the $|G(g_{i})\rangle$ is abbreviated as $|G_{i}\rangle$. We
have assumed that the control term $\hat{Y}$ can be expressed as a summation
over even operators of independent mode
$\hat{Y}=\sum_{\vec{k}}\hat{Y}_{\vec{k}}$. If $O(T)$ can be represented as a
sum of even operators over independent modes, such that
$\hat{O}(T)=\sum_{\vec{k}}\hat{O}_{\vec{k}}(T)$, then it follows that:

$\displaystyle\langle G_{i}|\hat{O}(T)\hat{Y}(t)|G_{i}\rangle$ (35)
$\displaystyle=$ $\displaystyle\sum_{\vec{k},\vec{k}^{\prime}}\langle
G_{i,\vec{k}}|\hat{O}_{\vec{k}}(T)\hat{Y}_{\vec{k}^{\prime}}(t)|G_{i,\vec{k}^{\prime}}\rangle$
$\displaystyle=$ $\displaystyle\sum_{\vec{k}^{\prime}\not=\vec{k}}\langle
G_{i,\vec{k}}|\hat{O}_{\vec{k}}(T)\hat{Y}_{\vec{k}^{\prime}}(t)|G_{i,\vec{k}^{\prime}}\rangle+\sum_{\vec{k}}\langle
G_{i,\vec{k}}|\hat{O}_{\vec{k}}(T)\hat{Y}_{\vec{k}}(t)|G_{i,\vec{k}}\rangle$
$\displaystyle=$ $\displaystyle\sum_{\vec{k},\vec{k}^{\prime}}\langle
G_{i,\vec{k}}|\hat{O}_{\vec{k}}(T)|G_{i,\vec{k}}\rangle\langle
G_{i,\vec{k}^{\prime}}|\hat{Y}_{\vec{k}^{\prime}}(t)|G_{i,\vec{k}^{\prime}}\rangle$
$\displaystyle-$ $\displaystyle\sum_{\vec{k}}\langle
G_{i,\vec{k}}|\hat{O}_{\vec{k}}(T)|G_{i,\vec{k}}\rangle\langle
G_{i,\vec{k}}|\hat{Y}_{\vec{k}}(t)|G_{i,\vec{k}}\rangle$ $\displaystyle+$
$\displaystyle\sum_{\vec{k}}\langle
G_{i,\vec{k}}|\hat{O}_{\vec{k}}(T)\hat{Y}_{\vec{k}}(t)|G_{i,\vec{k}}\rangle$
$\displaystyle=$ $\displaystyle\sum_{\vec{k}}\langle
G_{i,\vec{k}}|\hat{O}_{\vec{k}}(T)|\bar{G}_{i,\vec{k}}\rangle\langle\bar{G}_{i,\vec{k}}|\hat{Y}_{\vec{k}^{\prime}}(t)|G_{i,\vec{k}}\rangle$
$\displaystyle=$
$\displaystyle\sum_{\vec{k}}\langle\phi_{\vec{k}}(T)|\hat{O}_{1,\vec{k}}|\bar{\phi}_{\vec{k}}(T)\rangle\langle\bar{\phi}_{\vec{k}}(t)|\hat{Y}_{\vec{k}}|\phi_{\vec{k}}(t)\rangle.$

In deriving the penultimate line of Eq. (35), we have used the identity

$\displaystyle\mathbb{I}_{\vec{k}}=|G_{i,\vec{k}}\rangle\langle
G_{i,\vec{k}}|+|\bar{G}_{i,\vec{k}}\rangle\langle\bar{G}_{i,\vec{k}}|.$ (36)

Then, we can obtain Eq. (II.2). If $O(T)$ can be expressed as a product over
even operators of independent modes
$\hat{O}(T)=\prod_{\vec{k}}\hat{O}_{\vec{k}}(T)$, it can be calculated as

$\displaystyle\frac{\delta O_{2}(T)}{\delta g(t)}$ (37) $\displaystyle=$
$\displaystyle 2\imaginary\sum_{\vec{k}^{\prime}}\langle
G_{i,\vec{k}^{\prime}}|\hat{O}_{2,\vec{k}^{\prime}}(T)\hat{Y}_{\vec{k}^{\prime}}(t)|G_{i,\vec{k}^{\prime}}\rangle\prod_{\vec{k}\not=\vec{k}^{\prime}}\langle
G_{i,\vec{k}}|\hat{O}_{2,\vec{k}}(T)|G_{i,\vec{k}}\rangle$ $\displaystyle=$
$\displaystyle 2\imaginary\prod_{\vec{k}}\langle
G_{i,\vec{k}}|\hat{O}_{2,\vec{k}}(T)|G_{i,\vec{k}}\rangle\sum_{\vec{k}^{\prime}}\frac{\langle
G_{i,\vec{k}^{\prime}}|\hat{O}_{2,\vec{k}^{\prime}}(T)\hat{Y}_{\vec{k}^{\prime}}(t)|G_{i,\vec{k}^{\prime}}\rangle}{\langle
G_{i,\vec{k}^{\prime}}|\hat{O}_{2,\vec{k}^{\prime}}(T)|G_{i,\vec{k}^{\prime}}\rangle}$
$\displaystyle=$ $\displaystyle 2\imaginary\prod_{\vec{k}}\langle
G_{i,\vec{k}}|\hat{O}_{2,\vec{k}}(T)|G_{i,\vec{k}}\rangle$
$\displaystyle\times$ $\displaystyle\sum_{\vec{k}^{\prime}}\frac{\langle
G_{i,\vec{k}^{\prime}}|\hat{O}_{2,\vec{k}^{\prime}}(T)|\bar{G}_{i,\vec{k}^{\prime}}\rangle\langle\bar{G}_{i,\vec{k}^{\prime}}|\hat{Y}_{\vec{k}^{\prime}}(t)|G_{i,\vec{k}^{\prime}}\rangle}{\langle
G_{i,\vec{k}^{\prime}}|\hat{O}_{2,\vec{k}^{\prime}}(T)|G_{i,\vec{k}^{\prime}}\rangle}$
$\displaystyle=$ $\displaystyle
2\imaginary\prod_{\vec{k}}\langle\phi_{\vec{k}}(T)|\hat{O}_{2,\vec{k}}|\phi_{\vec{k}}(T)\rangle$
$\displaystyle\times$
$\displaystyle\sum_{\vec{k}^{\prime}}\frac{\langle\phi_{\vec{k}^{\prime}}(T)|\hat{O}_{2,\vec{k}^{\prime}}|\bar{\phi}_{\vec{k}^{\prime}}(T)\rangle\langle\bar{\phi}_{\vec{k}^{\prime}}(t)|\hat{Y}_{\vec{k}^{\prime}}|\phi_{\vec{k}^{\prime}}(t)\rangle}{\langle\phi_{\vec{k}^{\prime}}(T)|\hat{O}_{2,\vec{k}^{\prime}}|\phi_{\vec{k}^{\prime}}(T)\rangle}.$

We also use the identity in Eq.(36) in the above derivation.

##  Appendix B The numerical cost in neural network expansion

Figure 9:  (a) Computational time and (b) memory requirement for varying
neural network scales with spin lengths of $N=24$, $N=50$, and $N=100$. The
notation $n_{1}\times n_{2}$ denotes a network configuration with $n_{1}$ and
$n_{2}$ neurons in the first and second hidden layers, respectively. Each
point represents the average computational time or memory usage derived from
ten iterations of neuron additions to networks of identical scale. The dashed
lines indicate linear fits to the data points.

The parameter scale of expanding a DCNN can be quantified as $O((n_{1}+M\delta
n_{1})(n_{2}+M\delta n_{2}))$, where $M$ represents the number of expanding
the neural network as the network approaches convergence. For an intuitive
understanding, when the initial network size is substantially smaller than
$M$, the parameter growth simplifies to $O(M^{2})$. Conversely, for an
initially larger network relative to $M$, the complexity remains closer to
$O(n_{1}n_{2})$. For simplicity, we only consider the computational time and
the memory consumption of increasing the neurons one time, transitioning the
network from an $n_{1}\times n_{2}$ to an $(n_{1}+\delta
n_{1})\times(n_{2}+\delta n_{2})$ structure. where $n_{1}\times n_{2}$ denotes
a network with $n_{1}$ neurons in the first and $n_{2}$ in the second hidden
layers, respectively. The overall impact on computational time and memory can
be approximated by $\sum_{i}x_{i}y_{i}$, with $x_{i}$ representing the number
of iterations until successful neuron integration, and $y_{i}$ capturing the
associated time or memory costs. It’s noteworthy that computational duration
is subject to a variety of external factors, including CPU performance and
specific elements of the matrix. Assuming a constant spin length, the
variation in numerical costs linked to different network sizes is illustrated
in Fig. 9. With the parameter scale increases, we observe an almost linear
surge in computational time, indicative of the added complexity and data
processing requirements. Memory usage similarly rises, though its growth
trajectory is shallower compared to computational time, which suggests a
higher data throughput and computational complexity for larger networks. An
extension in spin length further accentuates these resource demands. However,
a quantifiable relationship between network size expansion and computational
resource consumption requires further investigation.

## References

  * Farhi  _et al._ [2001] E. Farhi, J. Goldstone, S. Gutmann, J. Lapan, A. Lundgren, and D. Preda, Science 292, 472 (2001).
  * Chen  _et al._ [2009] X. Chen, B. Zeng, Z.-C. Gu, B. Yoshida, and I. L. Chuang, Phys. Rev. Lett. 102, 220501 (2009).
  * Verstraete  _et al._ [2004] F. Verstraete, J. J. García-Ripoll, and J. I. Cirac, Phys. Rev. Lett. 93, 207204 (2004).
  * Sørensen  _et al._ [2010] A. S. Sørensen, E. Altman, M. Gullans, J. V. Porto, M. D. Lukin, and E. Demler, Phys. Rev. A 81, 061603 (2010).
  * Ljubotina  _et al._ [2022] M. Ljubotina, B. Roos, D. A. Abanin, and M. Serbyn, PRX Quantum 3, 030343 (2022).
  * Shettell and Markham [2020] N. Shettell and D. Markham, Phys. Rev. Lett. 124, 110502 (2020).
  * Frérot and Roscilde [2018] I. Frérot and T. Roscilde, Phys. Rev. Lett. 121, 020402 (2018).
  * Altherr and Yang [2021] A. Altherr and Y. Yang, Phys. Rev. Lett. 127, 060501 (2021).
  * Kessler  _et al._ [2014] E. M. Kessler, I. Lovchinsky, A. O. Sushkov, and M. D. Lukin, Phys. Rev. Lett. 112, 150802 (2014).
  * Hudelist  _et al._ [2014] F. Hudelist, J. Kong, C. Liu, J. Jing, Z. Ou, and W. Zhang, Nature Communications 5, 3049 (2014).
  * Marciniak  _et al._ [2022] C. D. Marciniak, T. Feldker, I. Pogorelov, R. Kaubruegger, D. V. Vasilyev, R. van Bijnen, P. Schindler, P. Zoller, R. Blatt, and T. Monz, Nature 603, 604 (2022).
  * Yin  _et al._ [2023] P. Yin, X. Zhao, Y. Yang, Y. Guo, W.-H. Zhang, G.-C. Li, Y.-J. Han, B.-H. Liu, J.-S. Xu, G. Chiribella, G. Chen, C.-F. Li, and G.-C. Guo, Nature Physics 19, 1122 (2023).
  * Tóth and Apellaniz [2014] G. Tóth and I. Apellaniz, Journal of Physics A: Mathematical and Theoretical 47, 424006 (2014).
  * Szigeti  _et al._ [2021] S. S. Szigeti, O. Hosten, and S. A. Haine, Applied Physics Letters 118, 140501 (2021).
  * Giovannetti  _et al._ [2004] V. Giovannetti, S. Lloyd, and L. Maccone, Science 306, 1330 (2004).
  * Giovannetti  _et al._ [2011] V. Giovannetti, S. Lloyd, and L. Maccone, Nature Photonics 5, 222 (2011).
  * Huang  _et al._ [2014] J. Huang, S. Wu, H. Zhong, and C. Lee, Annual Review of Cold Atoms and Molecules , 365 (2014).
  * Zhang  _et al._ [2019] C. Zhang, T. R. Bromley, Y.-F. Huang, H. Cao, W.-M. Lv, B.-H. Liu, C.-F. Li, G.-C. Guo, M. Cianciaruso, and G. Adesso, Phys. Rev. Lett. 123, 180504 (2019).
  * Giorda and Paris [2010] P. Giorda and M. G. A. Paris, Phys. Rev. Lett. 105, 020503 (2010).
  * Ollivier and Zurek [2001] H. Ollivier and W. H. Zurek, Phys. Rev. Lett. 88, 017901 (2001).
  * Note [1] If there are $N$ qubits in separable states, it is known that the uncertainty (that is, the inverse of the sensitivity) scales as ${\cal O}(N^{-1/2})$, which is called the standard quantum limit, while quantum physics allows the ultimate scalings being ${\cal O}(N^{-1})$, i.e., the Heisenberg scaling, in the absence of noise and ${\cal O}(N^{-3/4})$ in the presence of realistic decoherence.
  * Chu and Cai [2022] Y. Chu and J. Cai, Phys. Rev. Lett. 128, 200501 (2022).
  * Lücke  _et al._ [2011] B. Lücke, M. Scherer, J. Kruse, L. Pezzé, F. Deuretzbacher, P. Hyllus, O. Topic, J. Peise, W. Ertmer, J. Arlt, L. Santos, A. Smerzi, and C. Klempt, Science 334, 773 (2011).
  * Strobel  _et al._ [2014] H. Strobel, W. Muessel, D. Linnemann, T. Zibold, D. B. Hume, L. Pezzè, A. Smerzi, and M. K. Oberthaler, Science 345, 424 (2014).
  * Rahmani and Chamon [2011] A. Rahmani and C. Chamon, Phys. Rev. Lett. 107, 016402 (2011).
  * Zaletel  _et al._ [2021] M. P. Zaletel, A. Kaufman, D. M. Stamper-Kurn, and N. Y. Yao, Phys. Rev. Lett. 126, 103401 (2021).
  * Farooq  _et al._ [2015] U. Farooq, A. Bayat, S. Mancini, and S. Bose, Phys. Rev. B 91, 134303 (2015).
  * Reiter  _et al._ [2016] F. Reiter, D. Reeb, and A. S. Sørensen, Phys. Rev. Lett. 117, 040501 (2016).
  * Wei  _et al._ [2023] Z.-Y. Wei, D. Malz, and J. I. Cirac, Phys. Rev. Res. 5, L022037 (2023).
  * Motta  _et al._ [2020] M. Motta, C. Sun, A. T. K. Tan, M. J. O’Rourke, E. Ye, A. J. Minnich, F. G. S. L. Brandão, and G. K.-L. Chan, Nature Physics 16, 205 (2020).
  * Choi  _et al._ [2023] J. Choi, A. L. Shaw, I. S. Madjarov, X. Xie, R. Finkelstein, J. P. Covey, J. S. Cotler, D. K. Mark, H.-Y. Huang, A. Kale, H. Pichler, F. G. S. L. Brandão, S. Choi, and M. Endres, Nature 613, 468 (2023).
  * Léonard  _et al._ [2023] J. Léonard, S. Kim, J. Kwan, P. Segura, F. Grusdt, C. Repellin, N. Goldman, and M. Greiner, Nature 619, 495 (2023).
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Nature 551, 579 (2017).
  * Maslova  _et al._ [2019] N. S. Maslova, P. I. Arseyev, and V. N. Mantsevich, Scientific Reports 9, 3130 (2019).
  * Kato [1950] T. Kato, Journal of the Physical Society of Japan 5, 435 (1950).
  * Tabatabaei  _et al._ [2021] S. Tabatabaei, H. Haas, W. Rose, B. Yager, M. Piscitelli, P. Sahafi, A. Jordan, P. J. Poole, D. Dalacu, and R. Budakian, Phys. Rev. Appl. 15, 044043 (2021).
  * Xiang-Bin and Keiji [2001] W. Xiang-Bin and M. Keiji, Phys. Rev. Lett. 87, 097901 (2001).
  * Dziarmaga [2010] J. Dziarmaga, Advances in Physics 59, 1063 (2010).
  * Barankov and Polkovnikov [2008] R. Barankov and A. Polkovnikov, Phys. Rev. Lett. 101, 076801 (2008).
  * Dziarmaga [2005] J. Dziarmaga, Phys. Rev. Lett. 95, 245701 (2005).
  * Jozsa [1994] R. Jozsa, Journal of Modern Optics 41, 2315 (1994).
  * Nielsen and Chuang [2010] M. A. Nielsen and I. Chuang,  _Quantum computation and quantum information_ (Cambridge University Press, 2010).
  * del Campo  _et al._ [2012] A. del Campo, M. M. Rams, and W. H. Zurek, Phys. Rev. Lett. 109, 115703 (2012).
  * Sau and Sengupta [2014] J. D. Sau and K. Sengupta, Phys. Rev. B 90, 104306 (2014).
  * Damski [2014] B. Damski, Journal of Statistical Mechanics: Theory and Experiment 2014, P12019 (2014).
  * Deffner  _et al._ [2014] S. Deffner, C. Jarzynski, and A. del Campo, Phys. Rev. X 4, 021013 (2014).
  * Campbell  _et al._ [2015] S. Campbell, G. De Chiara, M. Paternostro, G. M. Palma, and R. Fazio, Phys. Rev. Lett. 114, 177206 (2015).
  * Guéry-Odelin  _et al._ [2019] D. Guéry-Odelin, A. Ruschhaupt, A. Kiely, E. Torrontegui, S. Martínez-Garaot, and J. G. Muga, Rev. Mod. Phys. 91, 045001 (2019).
  * Koch  _et al._ [2022] C. P. Koch, U. Boscain, T. Calarco, G. Dirr, S. Filipp, S. J. Glaser, R. Kosloff, S. Montangero, T. Schulte-Herbrüggen, D. Sugny, and F. K. Wilhelm, EPJ Quantum Technology 9, 19 (2022).
  * Doria  _et al._ [2011] P. Doria, T. Calarco, and S. Montangero, Phys. Rev. Lett. 106, 190501 (2011).
  * Hegerfeldt [2013] G. C. Hegerfeldt, Phys. Rev. Lett. 111, 260501 (2013).
  * Caneva  _et al._ [2009] T. Caneva, M. Murphy, T. Calarco, R. Fazio, S. Montangero, V. Giovannetti, and G. E. Santoro, Phys. Rev. Lett. 103, 240501 (2009).
  * Giannelli  _et al._ [2022] L. Giannelli, S. Sgroi, J. Brown, G. S. Paraoanu, M. Paternostro, E. Paladino, and G. Falci, Physics Letters A 434, 128054 (2022).
  * Giovannetti  _et al._ [2003] V. Giovannetti, S. Lloyd, and L. Maccone, Phys. Rev. A 67, 052109 (2003).
  * Levitin and Toffoli [2009] L. B. Levitin and T. Toffoli, Phys. Rev. Lett. 103, 160502 (2009).
  * Deffner and Campbell [2017] S. Deffner and S. Campbell, Journal of Physics A: Mathematical and Theoretical 50, 453001 (2017).
  * Bukov  _et al._ [2019] M. Bukov, D. Sels, and A. Polkovnikov, Phys. Rev. X 9, 011034 (2019).
  * Zhang  _et al._ [2023] M. Zhang, H.-M. Yu, and J. Liu, npj Quantum Information 9, 97 (2023).
  * Carleo and Troyer [2017] G. Carleo and M. Troyer, Science 355, 602 (2017).
  * Choo  _et al._ [2018] K. Choo, G. Carleo, N. Regnault, and T. Neupert, Phys. Rev. Lett. 121, 167204 (2018).
  * Choo  _et al._ [2020] K. Choo, A. Mezzacapo, and G. Carleo, Nature Communications 11, 2368 (2020).
  * Zhu  _et al._ [2022] Y. Zhu, Y.-D. Wu, G. Bai, D.-S. Wang, Y. Wang, and G. Chiribella, Nature Communications 13, 6222 (2022).
  * Carleo  _et al._ [2018] G. Carleo, Y. Nomura, and M. Imada, Nature Communications 9, 5322 (2018).
  * Day  _et al._ [2019] A. G. R. Day, M. Bukov, P. Weinberg, P. Mehta, and D. Sels, Phys. Rev. Lett. 122, 020601 (2019).
  * Wu  _et al._ [2015] N. Wu, A. Nanduri, and H. Rabitz, Phys. Rev. B 91, 041115 (2015).
  * Mi  _et al._ [2022a] X. Mi, M. Ippoliti, C. Quintana, A. Greene, Z. Chen, J. Gross, F. Arute, K. Arya, J. Atalaya, _et al._ , Nature 601, 531 (2022a).
  * Frey and Rachel [2022] P. Frey and S. Rachel, Science Advances 8, eabm7652 (2022).
  * Chen  _et al._ [2022] I.-C. Chen, B. Burdick, Y. Yao, P. P. Orth, and T. Iadecola, Phys. Rev. Res. 4, 043027 (2022).
  * Mi  _et al._ [2022b] X. Mi, M. Sonner, M. Y. Niu, K. W. Lee, B. Foxen, R. Acharya, I. Aleiner, T. I. Andersen, F. Arute, _et al._ , Science 378, 785 (2022b).
  * Turner  _et al._ [2018] C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Phys. Rev. B 98, 155134 (2018).
  * Wild  _et al._ [2021] D. S. Wild, D. Sels, H. Pichler, C. Zanoci, and M. D. Lukin, Phys. Rev. Lett. 127, 100504 (2021).
  * Bravyi  _et al._ [2020] S. Bravyi, A. Kliesch, R. Koenig, and E. Tang, Phys. Rev. Lett. 125, 260505 (2020).
  * Kairys  _et al._ [2020] P. Kairys, A. D. King, I. Ozfidan, K. Boothby, J. Raymond, A. Banerjee, and T. S. Humble, PRX Quantum 1, 020320 (2020).
  * Schiffer  _et al._ [2022] B. F. Schiffer, J. Tura, and J. I. Cirac, PRX Quantum 3, 020347 (2022).
  * Zhang  _et al._ [2021] D.-B. Zhang, G.-Q. Zhang, Z.-Y. Xue, S.-L. Zhu, and Z. D. Wang, Phys. Rev. Lett. 127, 020502 (2021).
  * Kraus [2011] B. Kraus, Phys. Rev. Lett. 107, 250503 (2011).
  * Schauss [2018] P. Schauss, Quantum Science and Technology 3, 023001 (2018).
  * Kim  _et al._ [2010] K. Kim, M.-S. Chang, S. Korenblit, R. Islam, E. E. Edwards, J. K. Freericks, G.-D. Lin, L.-M. Duan, and C. Monroe, Nature 465, 590 (2010).
  * Mostame and Schützhold [2008] S. Mostame and R. Schützhold, Phys. Rev. Lett. 101, 220501 (2008).
  * Cervera-Lierta [2018] A. Cervera-Lierta, Quantum 2, 114 (2018).
  * Kim  _et al._ [2023] Y. Kim, A. Eddins, S. Anand, K. X. Wei, E. van den Berg, S. Rosenblatt, H. Nayfeh, Y. Wu, M. Zaletel, K. Temme, and A. Kandala, Nature 618, 500 (2023).
  * Sachdev [1999] S. Sachdev, Physics World 12, 33 (1999).
  * Kitaev [2003] A. Kitaev, Annals of Physics 303, 2 (2003).
  * Wu [2012] N. Wu, Physics Letters A 376, 3530 (2012).
  * Brif  _et al._ [2010] C. Brif, R. Chakrabarti, and H. Rabitz, New Journal of Physics 12, 075008 (2010).
  * Note [2] Zheng Cheng, Meng-Jiao Lyu, Takayuki Myo, Hisashi Horiuchi, Hiroshi Toki, Zhong-Zhou Ren, Masahiro Isaka, Meng-Yun Mao, Hiroki Takemoto, Niu Wan, Wen-Long You, and Qing Zhao, to appear.
  * Rach  _et al._ [2015] N. Rach, M. M. Müller, T. Calarco, and S. Montangero, Phys. Rev. A 92, 062343 (2015).
  * Goetz  _et al._ [2016] R. E. Goetz, M. Merkel, A. Karamatskou, R. Santra, and C. P. Koch, Phys. Rev. A 94, 023420 (2016).
  * Chakrabarti  _et al._ [2008] R. Chakrabarti, R. Wu, and H. Rabitz, Phys. Rev. A 78, 033414 (2008).
  * Giagkiozis and Fleming [2015] I. Giagkiozis and P. Fleming, Information Sciences 293, 338 (2015).
  * Huang  _et al._ [2022a] H.-Y. Huang, R. Kueng, G. Torlai, V. V. Albert, and J. Preskill, Science 377, eabk3333 (2022a).
  * Białończyk and Damski [2020] M. Białończyk and B. Damski, Journal of Statistical Mechanics: Theory and Experiment 2020, 013108 (2020).
  * Wu [2020] N. Wu, Phys. Rev. E 101, 042108 (2020).
  * Cabrera and Jullien [1987] G. G. Cabrera and R. Jullien, Phys. Rev. B 35, 7062 (1987).
  * Roland and Cerf [2002] J. Roland and N. J. Cerf, Phys. Rev. A 65, 042308 (2002).
  * Wang  _et al._ [2020] Y. Wang, Y. Su, X. Chen, and C. Wu, Phys. Rev. Appl. 14, 044043 (2020).
  * Bera  _et al._ [2022] A. K. Bera, S. M. Yusuf, S. K. Saha, M. Kumar, D. Voneshen, Y. Skourski, and S. A. Zvyagin, Nature Communications 13, 6888 (2022).
  * Wang  _et al._ [2016] Y. Wang, J. Zhang, C. Wu, J. Q. You, and G. Romero, Phys. Rev. A 94, 012328 (2016).
  * Bason  _et al._ [2012] M. G. Bason, M. Viteau, N. Malossi, P. Huillery, E. Arimondo, D. Ciampini, R. Fazio, V. Giovannetti, R. Mannella, and O. Morsch, Nature Physics 8, 147 (2012).
  * Cimmarusti  _et al._ [2015] A. D. Cimmarusti, Z. Yan, B. D. Patterson, L. P. Corcos, L. A. Orozco, and S. Deffner, Phys. Rev. Lett. 114, 233602 (2015).
  * Li [2023] X. Li, Scientific Reports 13, 14734 (2023).
  * Vepsäläinen  _et al._ [2019] A. Vepsäläinen, S. Danilin, and G. S. Paraoanu, Science Advances 5, eaau5999 (2019).
  * Schrödinger [1935] E. Schrödinger, Naturwissenschaften 23, 807 (1935).
  * Wineland [2013] D. J. Wineland, Rev. Mod. Phys. 85, 1103 (2013).
  * Zurek [2001] W. H. Zurek, Nature 412, 712 (2001).
  * Ofek  _et al._ [2016] N. Ofek, A. Petrenko, R. Heeres, P. Reinhold, Z. Leghtas, B. Vlastakis, Y. Liu, L. Frunzio, S. M. Girvin, L. Jiang, M. Mirrahimi, M. H. Devoret, and R. J. Schoelkopf, Nature 536, 441 (2016).
  * Sun  _et al._ [2021] F.-X. Sun, S.-S. Zheng, Y. Xiao, Q. Gong, Q. He, and K. Xia, Phys. Rev. Lett. 127, 087203 (2021).
  * Dowling [2008] J. P. Dowling, Contemporary Physics 49, 125 (2008).
  * Huang  _et al._ [2022b] J. Huang, H. Huo, M. Zhuang, and C. Lee, Phys. Rev. A 105, 062456 (2022b).
  * Giovannetti  _et al._ [2006] V. Giovannetti, S. Lloyd, and L. Maccone, Phys. Rev. Lett. 96, 010401 (2006).
  * Lee [2006] C. Lee, Phys. Rev. Lett. 97, 150402 (2006).
  * Pezzé and Smerzi [2009] L. Pezzé and A. Smerzi, Phys. Rev. Lett. 102, 100401 (2009).
  * Martin  _et al._ [2013] M. J. Martin, M. Bishof, M. D. Swallows, X. Zhang, C. Benko, J. von Stecher, A. V. Gorshkov, A. M. Rey, and J. Ye, Science 341, 632 (2013).
  * Wasilewski  _et al._ [2010] W. Wasilewski, K. Jensen, H. Krauter, J. J. Renema, M. V. Balabas, and E. S. Polzik, Phys. Rev. Lett. 104, 133601 (2010).
  * Le Sage  _et al._ [2013] D. Le Sage, K. Arai, D. R. Glenn, S. J. DeVience, L. M. Pham, L. Rahn-Lee, M. D. Lukin, A. Yacoby, A. Komeili, and R. L. Walsworth, Nature 496, 486 (2013).
  * Leibfried  _et al._ [2005] D. Leibfried, E. Knill, S. Seidelin, J. Britton, R. B. Blakestad, J. Chiaverini, D. B. Hume, W. M. Itano, J. D. Jost, C. Langer, R. Ozeri, R. Reichle, and D. J. Wineland, Nature 438, 639 (2005).
  * Cesa and Pichler [2023] F. Cesa and H. Pichler, Phys. Rev. Lett. 131, 170601 (2023).
  * Xia  _et al._ [2023] J. Xia, X. Zhang, X. Liu, Y. Zhou, and M. Ezawa, Phys. Rev. Lett. 130, 106701 (2023).
  * Bombín  _et al._ [2023] H. Bombín, C. Dawson, R. V. Mishmash, N. Nickerson, F. Pastawski, and S. Roberts, PRX Quantum 4, 020303 (2023).
  * Mueller  _et al._ [2023] N. Mueller, J. A. Carolan, A. Connelly, Z. Davoudi, E. F. Dumitrescu, and K. Yeter-Aydeniz, PRX Quantum 4, 030323 (2023).
  * Salvia  _et al._ [2023] R. Salvia, M. Mehboudi, and M. Perarnau-Llobet, Phys. Rev. Lett. 130, 240803 (2023).
  * Zhou  _et al._ [2023] S. Zhou, S. Michalakis, and T. Gefen, PRX Quantum 4, 040305 (2023).
  * Faist  _et al._ [2023] P. Faist, M. P. Woods, V. V. Albert, J. M. Renes, J. Eisert, and J. Preskill, PRX Quantum 4, 040336 (2023).
  * Liu  _et al._ [2023] Q. Liu, Z. Hu, H. Yuan, and Y. Yang, Phys. Rev. Lett. 130, 070803 (2023).