# Work extractability from energy eigenstates under optimized local operations

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  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 examine the relationship between the second law of thermodynamics and the
energy eigenstates of quantum many-body systems that undergo cyclic unitary
evolution. Using a numerically optimized control protocol, we analyze how the
work extractability is affected by the integrability of the system. Our
findings reveal that, in nonintegrable systems the number of work-extractable
energy eigenstates converges to zero, even when the local control operations
are optimized. In contrast, in integrable systems, there are exponentially
many eigenstates from which positive work can be extracted, regardless of the
locality of the control operations. We numerically demonstrate that such a
strikingly different behavior can be attributed to the number of athermal
energy eigenstates. Our results provide insights into the foundations of the
second law of thermodynamics in isolated quantum many-body systems, which are
expected to contribute to the development of quantum many-body heat engines.

##  I Introduction

Figure 1:  Graphical description of the object of this study: We consider the
work extraction from energy eigenstates of the quantum many-body systems. We
investigate the dependence of the extraction on the locality of the control
operators and integrability of the system.

The second law of thermodynamics in isolated quantum many-body systems has
drawn significant attention in recent years, driven by a crucial question in
statistical physics: how does macroscopic irreversibility originate from
microscopically reversible dynamics [1, 2, 3, 4, 5, 6]? In the context of work
extraction, in particular, the exploration of the relationship between the
second law of thermodynamics and thermal pure quantum states, which are
indistinguishable from the Gibbs ensemble at a macroscopic level, constitutes
a fundamental problem in statistical mechanics [7, 8, 9, 10, 11, 12].

One of the most well-known expressions of the second law of thermodynamics is
embodied in Planck’s principle, which claims that it is impossible to extract
work from the Gibbs state via adiabatic cycles. Correspondingly, it is
understood to be unattainable to extract work from canonical ensemble through
any cyclic operations, which is called the passivity of the Gibbs ensemble
[13, 14, 15]. In contrast, there does not exist any no-go principle to
prevents work extraction from pure quantum states; the second law can be
violated if one is able to perform arbitrary unitary operations with arbitrary
precision. Note that this does not contradict the eigenstate thermalization
hypothesis [16, 17, 18, 1, 19, 20, 21, 22, 23, 24, 25], since it only states
the macroscopic indistinguishability from the thermal equilibrium.

Despite the great amount of effort by the existing works, our understanding on
the relationship between the second law of thermodynamics and pure quantum
many-body states is limited, in particular when all control operations are
subject to local constraint (See Fig. 1). A prior study by Ref. [7] has
investigated the connection between the work extractability from single energy
eigenstates that undergo a simple quench dynamics. Ref. [7] reveals that, in a
nonintegrable system, it is impossible with quench dynamics to extract work
from any energy eigenstate corresponding to a positive temperature, while we
can extract work from an exponentially large number of eigenstates if both the
initial and the quenched Hamiltonians are integrable. While this previous work
shows that the integrability of the model yields a qualitative difference in
the work extractability, the control operation is limited to a simple class of
quench dynamics. In order to address more general situations, one needs to
carefully examine the entire degrees of freedom in the control operations.

Control locality | Initial Hamiltonian $\hat{H}(0)$
---|---
Nonintegrable | Integrable
Local | - | $O(\exp(cL))$
Non-local | $O(\exp(cL))$ | $O(\exp(cL))$
Table 1:  The size scaling of $D_{\rm pos},$ the count of work-extractable
energy eigenstates within a fixed energy shell that corresponds to a positive
temperature. The count $D_{\rm pos}$ converges to zero under large system size
in nonintegrable systems if the control operation is local, while it grows
exponentially in other cases.

In the present work, we analyze the work extractability of energy eigenstates
subject to numerically optimized cyclic control operations. As we summarized
in Table 1, for nonintegrable systems, the number of work-extractable states
corresponding to positive temperatures converges to zero even with the
optimized protocol if the control operations are local, while we find that
there are exponentially many work-extractable states when control operations
are global. In sharp contrast, integrable systems allow work extraction from
exponentially many energy eigenstates, whether the control operations are
local or global. We further numerically demonstrate that such a difference in
work extractability can be attributed to the distribution of the entanglement
entropy (EE); integrable systems have exponentially many athermal states so
that we may reduce the energy of the state without increasing the EE, which is
a prohibited scenario in nonintegrable systems since almost all eigenstates
are expected to be thermal.

The remainder of the paper is organized as follows. In Sec. II, we describe
the setup and the algorithms to optimize the protocols extracting work. In
Sec. III, we present our main numerical results. Finally, we give the
conclusion and discussion in Sec. IV.

##  II Setup

###  II.1 Work extraction from energy eigenstates

We present the definition of work extractable eigenstates assuming cyclic
control operations. Let $\hat{H}$ be a Hamiltonian that satisfies
$\hat{H}|E_{\alpha}\rangle=E_{\alpha}|E_{\alpha}\rangle$ where the $\alpha$-th
eigenstate with energy $E_{\alpha}$ is denoted as $|E_{\alpha}\rangle$. Under
some unitary time evolution $U(t)$, we define the work extraction from the
$\alpha$-th eigenstate at time $t$ as

$\displaystyle W_{\alpha}(t):$ $\displaystyle=$
$\displaystyle\operatorname{Tr}\left[\hat{H}\rho_{\alpha}(0)\right]-\operatorname{Tr}\left[\hat{H}\rho_{\alpha}(t)\right]$
(1) $\displaystyle=$ $\displaystyle
E_{\alpha}-\operatorname{Tr}\left[\hat{H}U(t)|E_{\alpha}\rangle\langle
E_{\alpha}|U^{\dagger}(t)\right],$ (2)

where we have denoted the time evolved eigenstate as
$\rho_{\alpha}(t):=U(t)\ket{E_{\alpha}}\bra{E_{\alpha}}U^{\dagger}(t)$. In
particular, we assume that the unitary $U$ is discretized as

$\displaystyle U(t)=\prod_{n=1}^{N_{t}}\exp(-i\delta t\hat{H}(n\delta t)),$
(3)

where $\hat{H}(t)$ is a time-dependent Hamiltonian that controls the dynamics
of the system. Here, the time step is homogeneously given as $\delta t$, and
$N_{t}=t/\delta t$ is the number of discrete steps. In practice, we consider
an implementable set of local operations $\mathcal{B}=\\{\hat{O}_{i}\\}$ and
take their linear combination to constitute the time-dependent Hamiltonian as

$\displaystyle\hat{H}(t)=\sum_{\hat{O}_{i}\in\mathcal{B}}\gamma_{i}(t)\hat{O}_{i},$
(4)

where $\gamma_{i}(t)$ denotes the coefficient of the $i$-th operator at time
$t$.

We also introduce a metric to measure the work extractability. We indicate the
number of eigenstates in some energy shell from which positive work is
extracted as

$\displaystyle
D_{\mathrm{pos}}(t):=\left|\\{\ket{E_{\alpha}}~{}|~{}W_{\alpha}(t)\geq\varepsilon
L,~{}E_{\alpha}\in{\rm shell}\\}\right|,$ (5)

where the threshold $\varepsilon$ is introduced for the purpose of numerical
stability.

###  II.2 Optimization of control operations

In order to quantify the performance of control protocols, we introduce a
reward function. As an example of a smooth function that explicitly rewards
$W_{\alpha}>\varepsilon L$, we define the following:

$\displaystyle
r(t):=\sum_{E_{\alpha}\in\mathrm{shell}}\sigma_{a}(w_{\alpha}(t)-\varepsilon)+c(w_{\alpha}(t)-\delta)\theta(\delta-
w_{\alpha}(t)),$ (6)

where $w_{\alpha}(t):=W_{\alpha}(t)/L$ is the work density with $L$ being the
system size, $\sigma_{a}(x):=1/(1+\exp(-ax))$ is a sigmoid function,
$\theta(x)$ is the unit step function, and $c,\varepsilon,\delta$ are some
hyperparameters that determine the behavior of the reward function.
Concretely, $a$ in the sigmoid function controls the width of its ascending
segment, while $c$ regulates the slope of the linear function so that work
below a certain threshold $\delta$ is penalized. In the following, we fix the
hyperparameters as $a=30,c=0.1,\varepsilon=0.15,\delta=0.3$. While we expect
that the main result is not affected significantly by the choice of the reward
function, we leave this as a future work.

To reveal the work extractability under optimal control techniques, we
consider two optimization methodologies: the gradient-based algorithm and the
deep reinforcement learning (RL) algorithm. In the gradient-based algorithm we
compute the gradient $\partial r/\partial\gamma_{i}$ with a constraint such
that the Frobenius norm of the control Hamiltonian is bounded as
$\|\hat{H}(t)\|_{2}\leq C$, where the norm upper bound is fixed as
$C=\sqrt{2Ld}$ with $d=2^{L}$ being the full Hilbert space dimension. As an
implementable operation set, we consider $\mathcal{B}_{k}$ as a set of
translationally invariant and spatial inversion-symmetric operators that act
on at most $k$ contiguous sites. On the other hand, the RL algorithm
constructs a strategy to choose a unitary from a discrete set of unitaries
$\left\\{U_{m}\right\\}$ at each time step, so that we can maximize the
discounted reward expectation value that is estimated using a deep neural
network. In this work, $\left\\{U_{m}\right\\}$ is generated from
$\mathcal{B}_{2}$. See Appendix B for details of the optimization methods.

##  III Main Results

Figure 2:  The time evolution of $D_{\rm pos}$ generated by control protocols
optimized by gradient-based method, RL, and simple quench. The initial and
final Hamiltonians are (a) nonintegrable and (b) integrable. The results by
the gradient-based method with various operation set $\mathcal{B}_{k}$ are
shown by real lines, while the results from the RL-optimized protocol and
simple quench dynamics are shown by blue dashed and orange dotted lines,
respectively. Note that the RL algorithm performs optimization of the protocol
corresponding to $k=2$. The black dash-dotted lines represent the number of
eigenstates in the given energy shell. The system size is $L=12$. The
threshold is taken as $\varepsilon=0.15$, which is identical to the width of
the energy shell (see Appendix A for the dependence on the choice of
$\varepsilon$). The quench dynamics is performed under the Hamiltonian (7)
with $(h,g)=(0,1.5)$, which is also employed in Ref. [7].  Figure 3:  The
system-size scaling of $D_{\rm pos}$ generated by the control protocol
optimized by the gradient-based method. (a)–(d) The dynamics of
$D_{\mathrm{pos}}$ during the protocol for $L=10,12,14,16$. The integrability
of the initial/final Hamiltonians and the locality of the operation are (a)
nonintegrable and local ($k=4$), (b) integrable and local ($k=4$), (c)
nonintegrable and global ($k=L/2$), and (d) integrable and global ($k=L/2$).
The absence of the line for $L=16$ in panel (a) represents that there is no
energy eigenstate from which positive work is extracted. The scaling of
$D_{\rm pos}(t)$ at $t=1$ is summarized in (e). The inset of panel (e) shows
that the $D_{\rm pos}(t)$ increase exponentially with the system size except
for the nonintegrable system under local control.  Figure 4:  The relationship
between the work extractability and the EE. The vertical axis represents the
half-chain EE $S(\rho_{\alpha}(t))$, where the horizontal axis shows the
energy density. In panel (a), we show the results for eigenstates with
positive work extraction in the integrable case with local control ($k=4$). We
observe that the values of EE are all increased. Other panels (b)-(e) show
results for both positive and negative work extraction. Here each panel
represents (b) nonintegrable case under local control ($k=4$), (c) integrable
case under local control ($k=4$), (d) nonintegrable case under global control
($k=L/2$), and (e) integrable case under global control ($k=L/2$),
respectively. All data in this figure correspond to the data at $L=16$ and
$t=1$ in Fig. 3.  Figure 5:  The relationship between the EE and the extracted
work of the energy eigenstates by the optimized protocols. The vertical and
horizontal axes represent the changes of the EE and the extracted work
density, respectively. Each panel shows the result for (a) nonintegrable and
local, (b) integrable and local, (c) nonintegrable and global, and (d)
integrable and global control protocol. The red markers represent the
eigenstates counted as $D_{\mathrm{pos}}$. The histograms on the left of the
panels only show the distribution of the eigenstates corresponding to
$D_{\mathrm{pos}}$. All data in this figure correspond to the data at $L=16$
and $t=1$ in Fig. 3.

In this study, we adopt as the target Hamiltonian the one-dimensional quantum
Ising model under periodic boundary condition:

$\displaystyle\hat{H}=\sum_{l=1}^{L}\hat{\sigma}_{l}^{z}\hat{\sigma}_{l+1}^{z}+h\hat{\sigma}_{l}^{z}+g\hat{\sigma}_{l}^{x},$
(7)

where $\hat{\sigma}_{l}^{x,z}$ is the Pauli operator acting on the $l$-th
site, $h$ and $g$ are the strength of the longitudinal and transverse magnetic
fields, respectively. In the subsequent analysis, we exclusively consider
$(h,g)=(0.9045,0.809)$ as the nonintegrable case and $(h,g)=(0,0.5)$ as the
integrable case. Furthermore, we focus on the energy shell specified from
energy density $E/L\in[-0.25,-0.1]$, and we limit our discussion to the zero-
momentum sector and the inversion-symmetric sector.

###  III.1 System-size scaling of the optimized work extraction

Firstly, we show in Fig. 2 the time evolution of $D_{\rm pos}(t)$ under
various control protocols. It is consistent with Ref. [7] that the simple
quench dynamics in nonintegrable systems does not extract work at all.
Meanwhile, as we can see from Fig. 2(a), the optimization of the control
protocol allows us to extract work even in nonintegrable systems. It is
noteworthy that $D_{\rm pos}(t)$ increases as we increase the locality $k$ of
the control protocol. This is somewhat expected behavior, since with larger
$k$ we have higher expressibility in the time evolution unitary. As long as
the optimization method is successfully performed, we expect $D_{\rm pos}$ to
grow with $k$m which is indeed observed in our numerical results. In
particular, with $k=8$ we can extract work from the entire energy shell in
$L=12$.

While the gradient-based method and the RL method are different from each
other in a sense that the latter only performs discrete optimization, we find
that the complex procedure of the RL algorithm allows one to avoid the local
minima, so that the optimized protocol for $k=2$ achieves higher $D_{\rm pos}$
by the RL compared to the gradient-based method (see Fig. 2(b)). This is a
remarkable benefit of utilizing the RL algorithm, considering the fact that
the expressibility of the control unitary is limited to discrete operation
set. Note that it is also in agreement with previous works that the RL method
is expected to provide a powerful way to determine the quantum control for
many-body systems [26, 27]. However, here we focus on small or medium-size
systems to elaborate on the qualitative difference in the work extractability,
and therefore rather focus on the gradient-based method, since it is
numerically less demanding for such system sizes.

Next, we further investigate how the work extractability is affected by the
integrability of the Hamiltonian and the locality of control Hamiltonian. As
is shown in Fig. 3(a) and (b), the integrability strongly impacts the size
scaling when the control operations are local. In the nonintegrable case, we
find that $D_{\mathrm{pos}}$ becomes zero for $L=16$, and further expect that
this holds for larger system size $L$ as we further elaborate in the next
section. In contrast, under global control, $D_{\mathrm{pos}}$ increases along
with the system size $L$ regardless of the integrability, as shown in Fig.
3(c), (d). We summarize such behaviors in Fig. 3(e) by plotting the size
scaling of $D_{\mathrm{pos}}$ at $t=1$, which increases exponentially except
for the nonintegrable system under local operation.

We remark that the qualitative difference originating from the integrability
is exhibited not only under local control operations, but also under global
operations as well. Namely, we observe that the control time $t$ to achieve
some fixed value of $D_{\rm pos}$ remains constant in integrable systems,
while it takes longer in nonintegrable systems. It is an interesting open
question to seek how the required time $t$ scales with the system size.

###  III.2 Work extractability and athermal entanglement entropy

To understand the striking difference in the work extractability, we focus on
the number of athermal states, namely the energy eigenstates whose EE is
significantly lower than those of thermal pure quantum states. In
nonintegrable systems, it is known that almost all eigenstates are thermal,
i.e., the EE converges to that of the canonical ensemble in the thermodynamic
limit [28, 29]. Since the thermal EE increases monotonically with energy if
one focuses on the energy corresponding to a positive temperature, one must
reduce the EE to extract work from the system. However, it has been pointed
out that, in general it requires exponentially long time under local operation
to decrease the EE of the state [30, 31], and thus there is no efficient way
to extract work from any energy eigenstate. In contrast, such a property is
not generally present in integrable models; there can be exponentially many
eigenstates whose EE is lower than the value of thermal EE [32, 29]. This
means that there is no principle that prohibits one from extracting work from
the system even with local operations.

To examine our conjecture, we analyze the change in the EE before and after
the control operation as $\Delta S=S(\rho_{\alpha}(0))-S(\rho_{\alpha}(t))$
where $S$ is the half-chain EE. In Fig. 4(a), we illustrate that all work-
extractable energy eigenstates in integrable system encounter increase in the
EE under local control (Fig. 4(c) also shows data for non-work-extractable
states). In other words, we are allowed to increase the EE with the control
unitary. This does not necessarily require exponentially long control time,
and thus expected to be achievable.

We emphasize that such a scenario is not allowed in nonintegrable systems.
Figure 4(b) shows that, in nonintegrable systems, the EEs of most eigenstates
within the energy shell are distributed near the thermal EEs, and the
fluctuation is suppressed to be exponentially small 111In the analysis, we
calculate the thermal EE of pure states following Ref. [52]. The authors of
Ref. [52] calculated the EE of the random pure state with an energy
constraint, which they called _ergodic bipartition_ state.. This means that,
if one desires to extract work, one must reduce the EE of the state by either
employing global control, as shown in Fig. 4(d) and (e), or exponentially long
unitaries with local terms.

Figure 5 shows the distribution of work and entropy change before and after
the control operation of $t=1$. We can again see from Fig. 5(a) and (b) that
the EE increases under local operation, and therefore work can be extracted
from many eigenstates only in integrable systems. In contrast, Fig. 5(c) and
(d) reveal the EE can decrease under global operations, and therefore work
extraction can be realized even in nonintegrable systems.

We summarize that these findings are consistent with the proposed scenario
that, the qualitative difference in the scaling of athermal eigenstates is
directly related with the work extractability in integrable and nonintegrable
systems under local operations. Furthermore, these results also imply that we
expect $D_{\mathrm{pos}}$ to be zero even for $L>16$ in nonintegrable systems
under local control, because the number of athermal eigenstates converges to
zero. Meanwhile, this is not the case when one introduces demanding operations
such as global operations or exponentially long circuits; we can reduce the EE
so that it is possible to extract work.

##  IV Conclusion

In this work, we have utilized numerically optimized quantum control protocols
to analyze the work extractability from energy eigenstates of isolated quantum
many-body system. Under local control, we find that the integrability is
crucial to allow work extraction from exponentially many eigenstates.
Conversely, under global control, such a qualitative difference is not
observed when the evolution time is sufficiently long. By performing further
analysis on the EE, we further find a convincing argument that the work
extractability is related with the number of athermal eigenstates. Namely,
large fluctuation of the EE from the thermal EE in integrable systems is
crucial to allow work extraction, while such a mechanism is not present in
nonintegrable systems.

We envision two intriguing future directions. First, it is interesting to
explore what is the control time required for positive work extraction. As
pointed out in Sec. III.1, our fixed duration analysis indicate that positive
work extraction requires longer time under larger size in nonintegrable
systems. More detailed understanding on the scaling of the control time is an
important open problem. For instance, one may perform numerical investigation
for large-scale systems using variational methods such as tensor network [34,
35, 36] or artificial neural networks [37, 38, 39, 40, 41]. Another
interesting future direction is to study how general our findings hold among
various quantum many-body systems. We envision that the discrepancy between
the integral and nonintegral systems persists under more general local
Hamiltonians, which can be naturally expected especially for translationally
invariant models. While the current study has focused on an one-dimensional
model, it is intriguing to explore higher-dimensional systems.

##### Acknowledgments.—

The authors wish to thank fruitful discussion with Toshihiro Yada. 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 and JST Grant Number JPMJPF2221. T.S. is
supported by JSPS KAKENHI Grant Number JP19H05796, JST CREST Grant Number
JPMJCR20C1, Japan, and JST ERATO-FS Grant Number JPMJER2204, Japan. N.Y. and
T.S. are also supported by Institute of AI and Beyond of The University of
Tokyo. This work is supported by IBM Quantum. The RL is performed on AI
Bridging Cloud Infrastructure (ABCI) of National Institute of Advanced
Industrial Science and Technology (AIST).

##  Appendix A The threshold dependence of the system-size scaling

Figure A1:  The dependence of the threshold $\varepsilon$ in Eq. (5) on the
finite-size scaling of $D_{\mathrm{pos}}$. We fix the parameter in the reward
function (6) as the ones used in main part. Note that the absence of values
corresponding to $L=16$ in the panel (a) indicates that $D_{\mathrm{pos}}$ is
zero.

We discuss the relationship between the threshold $\varepsilon$ in Eq. (5) and
the finite-size scaling of $D_{\mathrm{pos}}$ that is defined in the main text
as

$\displaystyle
D_{\mathrm{pos}}:=\left|\\{\ket{E_{\alpha}}\in\mathrm{shell}|W_{\alpha}(t)\geq\varepsilon
L\\}\right|.$ (8)

Note that we fix the parameters in the reward function $a,c,\delta$ (6) as the
ones used in the main text.

In Fig. A1, we confirm that the main result as summarized in Table 1 in the
main text is robust under the variation of $\varepsilon.$ While we observe
that $D_{\rm pos}$ seems to saturate for $\varepsilon=0.175$, we argue that
this is an artifact due to the vanishing gradient in the reward function in
the gradient-based method, rather than the physical phenomena itself. Under
the current choice of the reward function, the gradient-based method is valid
when $\varepsilon$ does not exceed that of the energy shell width, and
therefore we have taken $\varepsilon=0.15$ in the main text.

##  Appendix B Optimization algorithms

Figure B2:  The abstract deep NN architecture used for the RL in the present work. At each time step $t_{i}$, the deep NN takes the time-step index $i$ as input. After intermediate computations by fully-connected layers, LSTM [42], and the dueling network [43], the deep NN outputs the estimate of the Q function, from which we determine the next action.  Linear + ReLU 1 | $600\rightarrow 1024$
---|---
Linear + ReLU 2 | $1024\rightarrow 1024$
LSTM | $(1024,7,1)\rightarrow 1024$
Dueling network | $1024\rightarrow 7$
Table 2:  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 one-hot representation of the time-step index. The input dimension $600$ corresponds to the number of time step. The output dimension $7$ corresponds to the dimension of the action space.  Reward discount $\eta$ |  | 0.997
---|---|---
Minibatch size |  | 608
Sequence length |  | 40
Optimizer |  | Adam [44]
Optimizer setting | learning rate | $10^{-4}$
| $\varepsilon$ | $10^{-3}$
| $\beta$ | $(0.9,0.999)$
Replay ratio |  | 1
Gradient norms clip |  | 80
Table 3:  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. [45] and its previous non-LSTM
version, Ref. [46], for the details of the hyperparameters. The other
parameters follow the ones in Ref. [47].

###  B.1 Gredient-based algorithm

We explain the details of the gradient-based algorithm employed in our study.
We will first discuss how to calculate the coefficients $\\{\gamma_{i}(t)\\}$
in the time-dependent Hamiltonian at each time step. Next, we introduce the
set of operators $\mathcal{B}=\\{\hat{O}_{i}\\}$ that constitutes the time-
dependent Hamiltonian.

####  B.1.1 Optimizing time-dependent Hamiltonian

Given the reward function defined as in Eq. (6) in the main text, we determine
the coefficients $\\{\gamma_{i}(t)\\}_{i}$ at each time step $t$ such that the
following optimization problem is solved:

$\displaystyle\begin{aligned} &\text{maximize}&\frac{dr}{dt}\\\ &\text{subject
to}&\|\hat{H}(t)\|_{2}\leq C.\end{aligned}$ (9)

Recall that in the gradient-vased algorithm, the time-evolving Hamiltonian
$\hat{H}(t)$ (4) in the main text consists of a linear combination of
Hermitian operators:

$\displaystyle\hat{H}(t)=\sum_{\hat{O}_{i}\in\mathcal{B}}\gamma_{i}(t)\hat{O}_{i},$
(10)

where $\mathcal{B}=\\{\hat{O}_{i}\\}_{i}$ denotes a set of Hermitian
operators. We optimize the coefficients $\gamma_{i}(t)$ to maximize $dr/dt$,
the derivative of the reward function (6), and obtain the optimized
Hamiltonian at each time step.

In the following, we derive the explicit representation of the
Karush–Kuhn–Tucker (KKT) conditions, the generalization of the Lagrange
conditions to the problem subjected to inequality constraints, for the problem
defined in Eqs. (9). First, consider the derivative of the reward function
(6):

$\displaystyle\frac{dr}{dt}$ $\displaystyle=\sum_{\alpha}\frac{\partial
r}{\partial w_{\alpha}}\frac{dw_{\alpha}}{dt}.$ (11)

The derivative of $w_{\alpha}$ in Eq. (11) is calculated as

$\displaystyle\begin{aligned}
\frac{dw_{\alpha}}{dt}=&L^{-1}\frac{d}{dt}\left[E_{\alpha}-\operatorname{Tr}\left[\hat{H}(\tau)\rho_{\alpha}(t)\right]\right]\\\
=&-L^{-1}\operatorname{Tr}\left[\hat{H}(\tau)\frac{dU(t)}{dt}\ket{E_{\alpha}}\bra{E_{\alpha}}U^{\dagger}(t)\right]\\\
&-L^{-1}\operatorname{Tr}\left[\hat{H}(\tau)U(t)\ket{E_{\alpha}}\bra{E_{\alpha}}\frac{dU^{\dagger}(t)}{dt}\right]\\\
=&iL^{-1}\operatorname{Tr}\left[\hat{H}(\tau)[\hat{H}(t),\rho_{\alpha}(t)]\right]\\\
=&iL^{-1}\sum_{i}\gamma_{i}(t)\operatorname{Tr}\left[\hat{H}(\tau)[\hat{O}_{i},\rho_{\alpha}(t)]\right]\\\
=&iL^{-1}\sum_{i}\gamma_{i}(t)\operatorname{Tr}\left[[\hat{H}(\tau),\hat{O}_{i}]\rho_{\alpha}(t)\right],\end{aligned}$
(12)

where $[A,B]=AB-BA$ is the commutator, and the substitution of the time-
dependent Hamiltonian (10) yields the fifth line. By substituting Eq. (12)
into Eq. (11), the derivative of the reward function reads

$\displaystyle\frac{dr}{dt}$ $\displaystyle=\sum_{i}\gamma_{i}(t)Y_{i}(t),$
(13)

where we have defined

$\displaystyle Y_{i}(t):=iL^{-1}\sum_{\alpha}\frac{\partial r}{\partial
w_{\alpha}}\operatorname{Tr}\left[[\hat{H}(\tau),\hat{O}_{i}]\rho_{\alpha}(t)\right].$
(14)

Note that the partial derivative of $r$ in Eq. (14) is calculated as

$\displaystyle\begin{aligned} \frac{\partial r}{\partial
w_{\alpha}}&=\frac{\partial}{\partial
w_{\alpha}}\left[\sigma_{a}(w_{\alpha}-\varepsilon)+c(w_{\alpha}-\delta)\theta(\delta-
w_{\alpha})\right]\\\
&=\begin{cases}a\sigma_{a}(w_{\alpha}-\varepsilon)\left\\{1-\sigma_{a}(w_{\alpha}-\varepsilon)\right\\}+c&(w_{\alpha}<\delta)\\\
a\sigma_{a}(w_{\alpha}-\varepsilon)\left\\{1-\sigma_{a}(w_{\alpha}-\varepsilon)\right\\}&(w_{\alpha}\geq\delta),\end{cases}\end{aligned}$
(15)

and Eq. (14) can be calculated using the expectation values regarding
$\\{\rho_{\alpha}(t)\\}_{\alpha}$ at each time step.

Substituting the time-dependent Hamiltonian (10) into the constraint in
problem (9) yields

$\displaystyle\begin{aligned} C^{2}&\geq\|\hat{H}(t)\|_{2}^{2}\\\
&=\operatorname{Tr}\left[\hat{H}^{\dagger}(t)\hat{H}(t)\right]\\\
&=\sum_{i,j}\gamma_{i}^{*}(t)\gamma_{j}(t)\operatorname{Tr}\left[\hat{O}_{i}^{\dagger}\hat{O}_{j}\right]\\\
&=Ld\sum_{i}|\gamma_{i}(t)|^{2},\end{aligned}$ (16)

where we assume that
$\operatorname{Tr}[\hat{O}_{i}^{\dagger}\hat{O}_{j}]=dL\delta_{ij}$. As we
later see in the next subsection, this assumption is valid for the set
$\mathcal{B}$ considered in this work.

We define the Lagrange function as

$\displaystyle\mathcal{L}:=-\sum_{i}\gamma_{i}(t)Y_{i}(t)+\lambda\left(Ld\sum_{i}|\gamma_{i}(t)|^{2}-C^{2}\right).$
(17)

The KKT condition for the optimization problem (9) is expressed as

$\displaystyle\frac{\partial\mathcal{L}}{\partial\gamma_{i}}=-Y_{i}(t)+2\lambda
Ld\gamma_{i}=0\ \text{for any $i$},$ (18a) $\displaystyle
Ld\sum_{i}|\gamma_{i}(t)|^{2}-C^{2}\leq 0,$ (18b)
$\displaystyle\lambda\left(Ld\sum_{i}|\gamma_{i}(t)|^{2}-C^{2}\right)=0,$
(18c) $\displaystyle\lambda\geq 0.$ (18d)

Regarding the constraint, we consider the two cases: (i)
$\|\hat{H}(t)\|_{2}=C$ and (ii) $\|\hat{H}(t)\|_{2}<C$. In the case (i), the
condition (18) becomes

$\displaystyle\gamma_{i}=\frac{Y_{i}(t)}{2\lambda Ld}\ \text{for any $i$},$
(19a) $\displaystyle Ld\sum_{i}|\gamma_{i}(t)|^{2}=C^{2},$ (19b)
$\displaystyle\lambda\geq 0.$ (19c)

Substituting Eq. (19a) into Eq. (19b), we obtain

$\displaystyle\lambda=\frac{\|Y(t)\|}{2C\sqrt{Ld}},$ (20)

where $\|Y(t)\|=\sqrt{\sum_{i}|Y_{i}(t)|^{2}}$. By combining Eq. (19a) with
Eq. (20), the coefficients $\gamma$ are expressed as

$\displaystyle\gamma_{i}(t)=\frac{CY_{i}(t)}{\sqrt{Ld}\|Y(t)\|}.$ (21)

At each time step, $\gamma_{i}(t)$ is calculated using the expectation values
of $\hat{H}(\tau)$ and $[\hat{H}(\tau),\hat{O}_{i}]$.

In the case (ii), $\lambda$ must be zero because of Eq. (18c). By substituting
$\lambda=0$ into Eq. (18a), we obtain the condition: $Y_{i}(t)=0$ for any $i$.
This equality holds when $\rho_{\alpha}(t)$ commutes with $\hat{H}(\tau)$ for
any $\alpha$, which results in the vanishing of the gradient of the reward
function (see Eq. (12)). Note that it is also natural to expect that the case
(ii) is essentially not relevant, since the optimization is done such that the
gradient is maximized.

In our analysis, the initial states are the energy eigenstates of
$\hat{H}(0)=\hat{H}(\tau)$. Therefore, $\rho_{\alpha}(t)$ commutes with
$\hat{H}(\tau)$ for any $\alpha$, leading to the vanishing of the gradient of
the reward function. To perform the gradient-based algorithm, we add a
perturbation to the dynamics at the beginning of the simulation. Specifically,
we perform a time evolution generated by $\sum_{l}\hat{\sigma}_{l}^{x}$ for a
small duration $\delta t=0.001$.

####  B.1.2 The operator set for gradient-based algorithm

In this section, we construct $\mathcal{B}_{k}$, which is a translationally
invariant, spatial inversion-symmetric and orthogonal operator set that act on
at most $k$ contiguous sites.

The first step is to construct a translationally invariant operator set
$\mathcal{A}_{k}$. Consider the $L$-qubit Pauli group $\mathcal{P}_{L}:=\\{\pm
1,\pm
i\\}\cdot\\{\hat{\sigma}^{0},\hat{\sigma}^{x},\hat{\sigma}^{y},\hat{\sigma}^{z}\\}^{\otimes
L}$ where $\hat{\sigma}^{a}~{}(a=0,x,y,z)$ are Pauli operators. For a given
Pauli operator $\hat{P}_{a}\in\mathcal{P}_{L}$ that acts on at most $k$
contiguous sites, we take a linear combination of translated operators such
that $\hat{Q}_{a}:=\sum_{l=0}^{L-1}\hat{T}^{l}\hat{P}_{a}\hat{T}^{-l}$, where
$\hat{T}$ is one-site translation operator, e.g.,
$T\hat{\sigma}_{l}^{x}T^{-1}=\hat{\sigma}_{l+1}^{x}$. Then, from the set of
operator $\\{\hat{Q}_{a}\\}_{a}$, we choose the elements of $\mathcal{A}_{k}$
so that there is no duplication. Note that we regard the operators as
identical such that only their global phases differ.

Next, we create an inversion-symmetric operator subset from $\mathcal{A}_{k}$,
where $\hat{R}$ is the spatial inversion operator. Namely, we remove
inversion-asymmetric elements from $\mathcal{A}_{k}$ and adopt
$\hat{Q}_{a}^{{}^{\prime}}=\left(\hat{Q}_{a}+\hat{R}\hat{Q}_{a}\hat{R}\right)/\sqrt{2}$
as the element of $\mathcal{B}_{k}$ in such a way that there is no redundancy.

Finally, we confirm that the elements in $\mathcal{B}_{k}$ satisfy the
orthogonality. Concretely, we straightforwardly obtain the following:

$\displaystyle\begin{aligned}
\operatorname{Tr}\left[\hat{Q}_{a}^{\dagger}\hat{Q}_{b}\right]&=\operatorname{Tr}\left[\sum_{l=0}^{L-1}\hat{T}^{l}\hat{P}_{a}^{\dagger}\hat{T}^{-l}\sum_{l^{\prime}=0}^{L-1}\hat{T}^{l^{\prime}}\hat{P}_{b}\hat{T}^{-l^{\prime}}\right]\\\
&=\sum_{l,l^{\prime}=0}^{L-1}\operatorname{Tr}\left[\hat{P}_{a}^{\dagger}\hat{T}^{-(l-l^{\prime})}\hat{P}_{b}\hat{T}^{l-l^{\prime}}\right]\\\
&=\sum_{l,l^{\prime}=0}^{L-1}d\delta_{a,b}=\begin{cases}Ld&(a=b)\\\ 0&(a\neq
b).\end{cases}\end{aligned}$ (22)

Here, we used the fact that we have chosen $\hat{Q}_{a}\in\mathcal{A}_{k}$
such that there is no redundancy, which implies that $\hat{P}_{a}$ does not
coincide with any $\hat{P}_{b\neq a}$ under any translation operation.
Following similar calculation, we also confirm that the inversion-symmetrized
elements $\hat{Q}_{a}^{{}^{\prime}}$ satisfy the orthogonality and the norm
$\|\hat{Q}_{a}^{{}^{\prime}}\|_{2}=\sqrt{Ld}$. As a result, we have verified
the orthogonality of $\mathcal{B}_{k}$ and the norm of elements
$\|\hat{Q}_{a}\|_{2}=\|\hat{Q}_{a}^{{}^{\prime}}\|_{2}=\sqrt{Ld}$.

###  B.2 Deep reinforcement learning

When adopting deep RL, we construct the protocol with a unitary sequence, in
which each element corresponds to the time evolution at each time step, chosen
from a fixed set of unitaries. The elements of the fixed set of unitaries are
generated by the following Hermitian operator set $\\{\hat{H}_{m}\\}_{m}$:

$\displaystyle\begin{aligned}
&\sum_{l}J\sigma^{z}_{l}\sigma^{z}_{l+1}+h_{\mathrm{I}}\sigma^{z}_{l},\
\sum_{l}J\sigma^{z}_{l}\sigma^{z}_{l+1}+h_{\mathrm{N}}\sigma^{z}_{l},\\\
&\sum_{l}\hat{\sigma}_{l}^{x}\hat{\sigma}_{l+1}^{y}+\hat{\sigma}_{l}^{y}\hat{\sigma}_{l+1}^{x},\
\sum_{l}\hat{\sigma}_{l}^{y}\hat{\sigma}_{l+1}^{z}+\hat{\sigma}_{l}^{z}\hat{\sigma}_{l+1}^{y},\\\
&\sum_{l}g_{\mathrm{I}}\sigma^{x}_{l},\ \sum_{l}g_{\mathrm{N}}\sigma^{x}_{l}\
\sum_{l}\hat{\sigma}_{l}^{y},\end{aligned}$ (23)

where $(J,\ h_{\mathrm{I}},\ h_{\mathrm{N}},\ g_{\mathrm{I}},\
g_{\mathrm{N}})=(1,\ 0,\ 0.9045,\ 0.5,\ 0.809)$. Note that the norm of these
terms satisfies the upper bound $\sqrt{2Ld}$ adopted in the main text. These
terms, which are linear combinations of elements in $\mathcal{B}_{2}$, are
also used in Ref. [26, 27]. In Sec. III.1, we set the time-step length as
$0.04$ and the number of time-step is $600$.

Deep reinforcement learning, specifically deep Q-learning, utilizes a NN to
approximate the following optimal action-value function Q [48, 49]:

$\displaystyle
Q^{*}(t_{i},m)=\max_{\pi}\mathbb{E}_{\pi}\left[r_{t_{i}}+\sum_{n=1}^{\infty}\eta^{n}r_{t_{i+n}}\middle|\hat{H}(t_{i})=\hat{H}_{m},\
\pi\right],$ (24)

which denotes the maximum sum of rewards $r_{t}$ discounted by $\eta$
($0<\eta<1$) in a stochastic policy that selects actions based on a
probability distribution as $\pi\left\lparen
m|t_{i}\right\rparen=\Pr\left\lparen m|t_{i}\right\rparen$. A potent variant
of Q-learning harnesses the capabilities of deep NNs to represent the action-
value function, which is hence referred to as deep RL algorithm [50, 51].

In this paper, we direct our focus toward a non-distributed implementation
[47] of a deep RL algorithm termed R2D2 [45]. R2D2 is a form of deep
Q-learning algorithm, and assumes that the agent can obtain partial
information about the state of the environment. Figure B2 shows the overall
picture. The details regarding the network structure are shown in Table 2 and
Table 3.

## References

  * Rigol  _et al._ [2008] M. Rigol, V. Dunjko, and M. Olshanii, Thermalization and its mechanism for generic isolated quantum systems, Nature 452, 854 (2008).
  * Iyoda  _et al._ [2017] E. Iyoda, K. Kaneko, and T. Sagawa, Fluctuation theorem for Many-Body pure quantum states, Physical review letters 119, 100601 (2017).
  * Tasaki [2000] H. Tasaki, The second law of thermodynamics as a theorem in quantum mechanics, arXiv:cond-mat/0011321 [cond-mat.stat-mech] (2000).
  * Goldstein  _et al._ [2013] S. Goldstein, T. Hara, and H. Tasaki, The second law of thermodynamics for pure quantum states, arXiv:1303.6393 [cond-mat.stat-mech] (2013).
  * Ikeda  _et al._ [2015] T. N. Ikeda, N. Sakumichi, A. Polkovnikov, and M. Ueda, The second law of thermodynamics under unitary evolution and external operations, Annals of physics 354, 338 (2015).
  * Kaneko  _et al._ [2017] K. Kaneko, E. Iyoda, and T. Sagawa, Saturation of entropy production in quantum many-body systems, Physical review. E 96, 062148 (2017).
  * Kaneko  _et al._ [2019] K. Kaneko, E. Iyoda, and T. Sagawa, Work extraction from a single energy eigenstate, Physical review. E 99, 032128 (2019).
  * Dorner  _et al._ [2013] R. Dorner, S. R. Clark, L. Heaney, R. Fazio, J. Goold, and V. Vedral, Extracting quantum work statistics and fluctuation theorems by single-qubit interferometry, Physical review letters 110, 230601 (2013).
  * Gallego  _et al._ [2014] R. Gallego, A. Riera, and J. Eisert, Thermal machines beyond the weak coupling regime, New journal of physics 16, 125009 (2014).
  * Perarnau-Llobet  _et al._ [2016] M. Perarnau-Llobet, A. Riera, R. Gallego, H. Wilming, and J. Eisert, Work and entropy production in generalised gibbs ensembles, New journal of physics 18, 123035 (2016).
  * Modak and Rigol [2017] R. Modak and M. Rigol, Work extraction in an isolated quantum lattice system: Grand canonical and generalized gibbs ensemble predictions, Physical review. E 95, 062145 (2017).
  * Le  _et al._ [2018] T. P. Le, J. Levinsen, K. Modi, M. M. Parish, and F. A. Pollock, Spin-chain model of a many-body quantum battery, Physical review. A 97, 022106 (2018).
  * Pusz and Woronowicz [1978] W. Pusz and S. L. Woronowicz, Passive states and KMS states for general quantum systems, Communications in Mathematical PHysics 58, 273 (1978).
  * Lenard [1978] A. Lenard, Thermodynamical proof of the gibbs formula for elementary quantum systems, Journal of statistical physics 19, 575 (1978).
  * Mitsuhashi  _et al._ [2022] Y. Mitsuhashi, K. Kaneko, and T. Sagawa, Characterizing symmetry-protected thermal equilibrium by work extraction, Phys. Rev. X 12, 021013 (2022).
  * Neumann [1929] J. v. Neumann, Beweis des ergodensatzes und desH-Theorems in der neuen mechanik, Zeitschrift für Physik 57, 30 (1929).
  * Deutsch [1991] J. M. Deutsch, Quantum statistical mechanics in a closed system, Physical review. A 43, 2046 (1991).
  * Srednicki [1994] M. Srednicki, Chaos and quantum thermalization, Physical review. E, Statistical physics, plasmas, fluids, and related interdisciplinary topics 50, 888 (1994).
  * 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).
  * Beugeling  _et al._ [2014] W. Beugeling, R. Moessner, and M. Haque, Finite-size scaling of eigenstate thermalization, Physical Review E 89, 042112 (2014).
  * Yoshizawa  _et al._ [2018] T. Yoshizawa, E. Iyoda, and T. Sagawa, Numerical Large Deviation Analysis of the Eigenstate Thermalization Hypothesis, Physical Review Letters 120, 200604 (2018).
  * Steinigeweg  _et al._ [2013] R. Steinigeweg, J. Herbrych, and P. Prelovšek, Eigenstate thermalization within isolated spin-chain systems, Physical Review E 87, 012118 (2013).
  * Sorg  _et al._ [2014] S. Sorg, L. Vidmar, L. Pollet, and F. Heidrich-Meisner, Relaxation and thermalization in the one-dimensional Bose-Hubbard model: A case study for the interaction quantum quench from the atomic limit, Physical Review A 90, 033606 (2014).
  * Khodja  _et al._ [2015] A. Khodja, R. Steinigeweg, and J. Gemmer, Relevance of the eigenstate thermalization hypothesis for thermal relaxation, Physical Review E 91, 012120 (2015).
  * Mondaini  _et al._ [2016] R. Mondaini, K. R. Fratus, M. Srednicki, and M. Rigol, Eigenstate thermalization in the two-dimensional transverse field Ising model, Physical Review E 93, 032104 (2016).
  * 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).
  * Baba  _et al._ [2023] S. Z. Baba, N. Yoshioka, Y. Ashida, and T. Sagawa, Deep reinforcement learning for preparation of thermal and prethermal quantum states, Physical Review Applied , 014068 (2023).
  * Deutsch  _et al._ [2013] J. M. Deutsch, H. Li, and A. Sharma, Microscopic origin of thermodynamic entropy in isolated systems, Physical review. E, Statistical, nonlinear, and soft matter physics 87, 042135 (2013).
  * Beugeling  _et al._ [2015] W. Beugeling, A. Andreanov, and M. Haque, Global characteristics of all eigenstates of local many-body hamiltonians: participation ratio and entanglement entropy, Journal of statistical mechanics 2015, P02002 (2015).
  * Metz and Bukov [2022] F. Metz and M. Bukov, Self-Correcting quantum Many-Body control using reinforcement learning with tensor networks, arXiv:2201.11790 [quant-ph] (2022).
  * Poulin  _et al._ [2011] D. Poulin, A. Qarry, R. Somma, and F. Verstraete, Quantum simulation of time-dependent hamiltonians and the convenient illusion of hilbert space, Physical review letters 106, 170501 (2011).
  * Alba  _et al._ [2009] V. Alba, M. Fagotti, and P. Calabrese, Entanglement entropy of excited states, Journal of statistical mechanics 2009, P10020 (2009).
  * Note [1] In the analysis, we calculate the thermal EE of pure states following Ref. [52]. The authors of Ref. [52] calculated the EE of the random pure state with an energy constraint, which they called _ergodic bipartition_ state.
  * Vidal [2003] G. Vidal, Efficient classical simulation of slightly entangled quantum computations, Phys. Rev. Lett. 91, 147902 (2003).
  * Vidal [2004] G. Vidal, Efficient simulation of one-dimensional quantum many-body systems, Phys. Rev. Lett. 93, 040502 (2004).
  * Daley  _et al._ [2004] A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal, Time-dependent density-matrix renormalization-group using adaptive effective hilbert spaces, Journal of Statistical Mechanics: Theory and Experiment 2004, P04005 (2004).
  * Carleo and Troyer [2017] G. Carleo and M. Troyer, Solving the quantum many-body problem with artificial neural networks, Science 355, 602 (2017).
  * Yoshioka and Hamazaki [2019] N. Yoshioka and R. Hamazaki, Constructing neural stationary states for open quantum many-body systems, Phys. Rev. B 99, 214306 (2019).
  * Hartmann and Carleo [2019] M. J. Hartmann and G. Carleo, Neural-network approach to dissipative quantum many-body dynamics, Phys. Rev. Lett. 122, 250502 (2019).
  * 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).
  * Vicentini  _et al._ [2019] F. Vicentini, A. Biella, N. Regnault, and C. Ciuti, Variational neural-network ansatz for steady states in open quantum systems, Phys. Rev. Lett. 122, 250503 (2019).
  * Hochreiter and Schmidhuber [1997] S. Hochreiter and J. Schmidhuber, Long Short-Term Memory, Neural Computation 9, 1735 (1997).
  * 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.
  * Kingma and Ba [2015] D. P. Kingma and J. Ba, Adam: A Method for Stochastic Optimization, in  _International Conference on Learning Representations_ (2015).
  * 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).
  * 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).
  * [47] A. Stooke and P. Abbeel, rlpyt: A Research Code Base for Deep Reinforcement Learning in PyTorch, arXiv:1909.01500 .
  * 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).
  * 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).
  * Li [2018] Y. Li, Deep reinforcement learning, arXiv:1810.06339 [cs.LG] (2018).
  * Henderson  _et al._ [2017] P. Henderson, R. Islam, P. Bachman, J. Pineau, D. Precup, and D. Meger, Deep reinforcement learning that matters, arXiv:1709.06560 [cs.LG] (2017).
  * Lu and Grover [2019] T.-C. Lu and T. Grover, Renyi entropy of chaotic eigenstates, Physical review. E 99, 032111 (2019).


# Leveraging Staggered Tessellation for Enhanced Spatial Resolution in High-
Granularity Calorimeters

Sebouh J. Paul  Miguel Arratia

###### Abstract

We advance the concept of high-granularity calorimeters with staggered
tessellations, underscoring the effectiveness of a design incorporating
multifold staggering cycles based on hexagonal cells to enhance position
resolution. Moreover, we introduce HEXPLIT, a sub-cell re-weighting algorithm
tailored to harness staggered designs, resulting in additional performance
improvements. By combining our proposed staggered design with HEXPLIT, we
achieve an approximately twofold enhancement in position resolution for
neutrons across a wide energy range, as compared to unstaggered designs. These
findings hold the potential to elevate particle-flow performance across
various forthcoming facilities.

##  1 Introduction

Highly-granular calorimeters, optimized for particle-flow applications [1, 2],
offer unprecedented potential for the detailed measurement of particle
showers. Hits in multiple layers can be used to determine the axis of a
particle shower or the trajectory of a minimum-ionizing particle (MIP) with
great accuracy.

When a MIP hits two cells with the same transverse position, the second hit
adds very little information about the transverse position of its trajectory.
However, if the cells are staggered and only partially overlap with each
other, it becomes more likely that the particle passed through the overlapping
region. Exploiting this fact in the calorimeter’s design by staggering cell
positions in alternating layers improves its spatial resolution. The same
principle applies to particle showers, provided that the cell size is
comparable to the calorimeter’s Moliere radius.

Staggering detector elements is a common technique in calorimeter design aimed
at reducing dead areas and minimizing channeling effects. For instance, in the
CMS’s HGCAL design [3], the endcap petals are staggered with azimuthal
rotations within three groups of layers to eliminate dead areas between petals
[4]. Additionally, in the ATLAS TileCal design [5], quadrangular scintillator
cells are positioned perpendicular to the beamline and staggered in the
longitudinal direction.

In specific instances, staggering has been deliberately incorporated into
detector designs to enhance spatial resolution. For example, the CALICE
silicon-tungsten ECAL [6, 7] employs staggering of square cells to improve
spatial resolution, as shown by a test beam study [8]. On the other hand, the
CALICE tungsten-scintillator ECAL [9, 10] employs a configuration of
alternating vertical and horizontal strips.

While staggering alone can enhance the performance of traditional shower
reconstruction algorithms, it also opens up possibilities for specialized
algorithms that exploit the cell overlap. For instance, the Strip-Splitting
Algorithm (SSA) [11] was developed for calorimeters featuring alternating
vertical and horizontal strips. This algorithm partitions each strip into sub-
cells, to which energy is allocated based on the energy distribution in
adjacent layers. As a result, the performance achieved with strips closely
approximates what would be attainable with square pixels of dimensions
equivalent to the strip widths.

In this study, we advance the concept of staggering layers within high
granularity calorimeters to enhance spatial resolution. Specifically, we
introduce a novel staggering design based on a multifold repeating
tessellation patterns. Additionally, we introduce the HEXPLIT algorithm, which
harnesses the overlap between cells, akin to the way SSA utilizes strip
layouts. What sets HEXPLIT apart from SSA is its use of overlaps between
multiple layers, including adjacent and next-to-neighboring layers, each
featuring distinct layouts. As a result, HEXPLIT can capitalize on multifold
repeating patterns, yielding elevated levels of “effective granularity”.

This paper is organized as follows: we present multifold repeating
tessellation patterns in Section 2, outline our simulation framework in
Section 3, introduce the HEXSPLIT algorithm in Section 4, discuss our findings
in Section 5, and conclude in Section 6.

##  2 Tessellations

The choice of cell shape for creating a tessellation of an area influences the
detector performance and potential enhancements of staggering configurations.
The area, $A$, of a cell determines the total number of channels and the
detector granularity. The circumradius, $r_{c}$, or maximum distance between
the center of the cell and the vertices, is relevant for the light yield and
cell uniformity when using the SiPM-on-tile technology [12, 13]. The amount of
dead area between cells is proportional to the perimeter, $p$, of the cells.
The single-cell spatial resolution is related to the RMS displacement from the
centroid of the cell, defined geometrically by:

$\Delta x=\sqrt{\frac{\iint_{A}dx\,dy\,x^{2}}{\iint_{A}dx\,dy}};\hskip
28.45274pt\Delta y=\sqrt{\frac{\iint_{A}dx\,dy\,y^{2}}{\iint_{A}dx\,dy}}.$
(2.1)

For squares and rectangular strips, this yields the well-known result of
$1/\sqrt{12}$ times the width or length. It should be noted that the
resolution for showering particles can surpass the value given in Eq. 2.1 due
to energy being distributed across cells. Nevertheless, this remains a
valuable benchmark.

Table 1 presents a summary of various properties pertaining to tessellating
cell shapes that hold significance in calorimetry. It quantifies the well-
known fact that hexagons possess the smallest perimeter and the smallest
circumradius among all tessellating polygons for a given area. Moreover,
hexagons have the highest count of symmetry axes, thereby enabling intricate
staggering patterns.

| reg. triangle | square | reg. hexagon
---|---|---|---
$A$ | $\frac{\sqrt{3}}{4}s^{2}\approx 0.433s^{2}$ | $s^{2}$ | $\frac{3\sqrt{3}}{2}s^{2}\approx 2.60s^{2}$
$r_{c}$ | $\frac{2}{3^{3/4}}\sqrt{A}\approx 0.877\sqrt{A}$ | $\frac{1}{\sqrt{2}}\sqrt{A}\approx 0.707\sqrt{A}$ | $\sqrt{\frac{2}{3\sqrt{3}}}\sqrt{A}\approx 0.620\sqrt{A}$
$p$ | $\frac{6}{\sqrt{\sqrt{3}}}\sqrt{A}\approx 4.46\sqrt{A}$ | $4\sqrt{A}$ | $6\sqrt{\frac{2}{3\sqrt{3}}}\sqrt{A}\approx 3.72\sqrt{A}$
$\Delta x$ | $\frac{1}{\sqrt{24}}s\approx 0.204s$ | $\frac{1}{\sqrt{12}}s\approx 0.289s$ | $\sqrt{\frac{5}{24}}s\approx 0.456s$
$\approx 0.310\sqrt{A}$ | $\approx 0.289\sqrt{A}$ | $\approx 0.283\sqrt{A}$
$\Delta y$ | same as $\Delta x$ | same as $\Delta x$ | same as $\Delta x$
Table 1: Summary of the properties of tessellating regular polygons, including
the area ($A$), side length ($s$), circumradius ($r_{c}$), perimeter ($p$),
and the RMS displacement in the $x$ ($\Delta x$) and $y$ ($\Delta y$)
directions.

Figure 1 displays various tessellations, both square- and hexagon-based, with
and without staggering. The top row illustrates a square tessellation without
staggering on the left and with a simple111In the CALICE ECAL beam test, a
more intricate staggering pattern was employed along one axis [6]. staggering
pattern on the right, denoted as S2. Additionally, we introduce two novel
staggered tessellation patterns with hexagonal cells, labeled as H3 and H4.
These two staggering configurations repeat every three and four layers,
respectively.

Figure 1: Top row, from left to right: unstaggered squares, and the S2
iteration of the square tessellation. Bottom row, from left to right: non-
staggered hexagons, the H3 iteration of the hexagonal tessellation, and the H4
iteration of the hexagonal tessellation.

In H3, the centers of the cells on one layer have the same transverse position
as the corners of the cells on the other layers. Equivalently, the centers of
the cells in subsequent layers are one hexagon radius apart. This produces
triangular subcells, each of which is 1/6 of the area of the hexagonal cells
and has $\Delta x_{\rm tri}=\Delta y_{\rm tri}\approx 0.204s_{\rm hex}$, where
$s_{\rm hex}$ is the side length of the hexagons.

In H4, the centers of the cells on one layer are at the midpoints of the edges
of the cells in the other layers. Equivalently, the cell centers are shifted
by $\sqrt{3}/2$ hexagon radii from one layer to another. Effectively, groups
of four layers define rhombus-shaped subcells, each of which is 1/12 of the
area of the hexagonal cells. For a rhombus, the RMS displacement from the
centroid of the cell is $\frac{1}{\sqrt{8}}s\approx 0.354s$ in the long
direction and $\frac{1}{\sqrt{24}}s\approx 0.204s$ in the short direction,
where $s$ is the length of one side of the rhombus. However, since the side
length of the rhombical sub-cells is half of that of a hexagonal cell, we have
$\Delta x_{\rm rhombus}\approx 0.177s_{\rm hex}$ and $\Delta y_{\rm
rhombus}\approx 0.102s_{\rm hex}$.

All of these tessellation patterns could be practically implemented in
realistic detector designs. For example, they could be seamlessly integrated
into the SiPM-on-tile approach, where extensive layers of scintillator cells
are tessellated to span a large-area surface. This integration could be
achieved by encasing individual cells in reflective foil, akin to the CALICE
AHCAL [14] and CMS HGCAL [3] configurations. Additional techniques, like
attaching individual cells onto plastic frames [15] to establish unique
patterns, could also present viable alternatives.

##  3 Simulations

We employed the DD4HEP [16] framework, which internally relies on GEANT4 [17]
with the FTFP_BERT_HP physics list, to simulate the geometry222We adopted a
coordinate system aligned with the calorimeter, wherein the $z$-axis is
perpendicular to the layers, $y$ represents the vertical direction, and $x$ is
orthogonal to both $y$ and $z$, establishing a right-handed coordinate system.
of an iron-scintillator calorimeter akin to the CALICE AHCAL [14], CMS HGCAL
[3], and the ePIC forward calorimeter insert [18]. This setup comprises
alternating layers of 20 mm-thick iron absorbers and layers of 3 mm-thick
scintillator tiles. In total, 64 layers were included, for a total of
8.0$\lambda$. The detector area is 6$\times$6 m2.

Each cell is represented as either a hexagonal or square shape, both with an
area of 25 cm2. The circumradius of the cells is approximately 1.0 or 1.1
times the Molière radius of this detector (3.1 cm), respectively for the
hexagonal and square shapes. A repeating tessellation pattern using either
hexagons or square cells is implemented. Since DD4HEP does not include built-
in code for hexagonal segmentation, we developed our own plugin for this
purpose [19].

We simulated 2,000 single-neutron events within the energy range of 10–300 GeV
for each configuration (unstaggered squares, S2, unstaggered hexagons, H3, and
H4), covering polar angles from 0 to 5.5 mrad to avoid potential artifacts
that might arise from constraining the angle to a single value. The part of
the detector’s front face impacted by neutrons encompassed an area of
approximately 1200 cm2. Furthermore, we conducted additional simulations
involving muons and electrons to determine the hit-energy MIP unit and
establish the calorimeter’s electromagnetic energy scale. The MIP unit,
defined as the most probable value, was found to be MIP$=0.5$ MeV.

Hits are reconstructed using a 15-bit ADC that is free from electronic noise,
as well as optical and electronic cross-talk. In the subsequent analysis, only
hits with $E>0.1$ MIP and $t<200$ ns are taken into account. All hits
surpassing these thresholds are retained for further analysis, without
applying any clustering algorithms.

##  4 Shower-Reconstruction Algorithm

As a baseline, we evaluate the single-particle position resolution using a
conventional procedure, as described in Section 4.1. Subsequently, we
introduce HEXPLIT, an algorithm that exploits sub-cell information, outlined
in Section 4.2.

###  4.1 Log-weighted determination of the shower position

The position of the shower is reconstructed using a logarithmically weighted
average of the center positions of the cells 333We also experimented with both
the baseline and HEXPLIT algorithms by employing weights linearly proportional
to energy in neutron showers. However, this approach yielded significantly
poorer resolutions compared to the logarithmic weighting discussed here.,
similar to Ref. [20, 21]. We used the following formula:

$\vec{x}_{\rm recon}=\frac{\sum\limits_{i\in\rm
hits}\vec{x}_{i}w_{i}}{\sum\limits_{i\in\rm hits}w_{i}}$ (4.1)

where $\vec{x}_{i}$ are the 3d positions of the hits, and the weights,
$w_{i}$, are determined by

$w_{i}=\max\left(0,w_{0}+\ln{\frac{E_{i}}{E_{\rm tot}}}\right)$ (4.2)

Here, $E_{i}$ represents the energy of an individual hit, and $E_{\rm tot}$
denotes the total energy of the shower. The parameter $w_{0}$ governs the
minimum hit energy used. The optimal value for $w_{0}$ varies depending on
energy and hardware parameters, such as cell size and the Molière radius of
the detector.

###  4.2 Sub-cell reweighting with HEXPLIT

The HEXPLIT algorithm involves the energy reweighting of sub-cells formed
through the overlap of cells within a specific layer and neighboring layers.
The initial step involves characterizing the energy distribution within these
sub-cells. Following that, this process uses the positions and energies of
these sub-cells as input for the logarithmic weighting procedure outlined in
Section 4.1.

Assigning a weight to each sub-cell in a staggering cycle, denoted by the
number of layouts $N$, is accomplished using the following formula:

$W_{i}=\prod\limits_{j=1}^{N-1}{\rm max}(E_{j},\delta),$ (4.3)

Here, $i$ denotes the index of a subcell. The product over $j$ is over the
tiles in the nearby layers that overlap with subcell $i$, and $E_{j}$
represents the energy in that tile. The parameter $\delta$ acts as a noise
threshold and prevents division by zero. For this analysis, we set $\delta$ to
1 MIP.

In scenarios where $N$ is even, such as in the instances of H4 and S2, the
layouts on layers that are $N/2$ layers distant from the current layer exhibit
identical patterns. Consequently, we introduce a virtual tile, whose energy is
the aggregate of the energies from the tiles within those layers overlapping
the $i^{\rm th}$ subcell.

In the H3 configuration, this is analogous to taking a product over the
overlapping tiles in the layer immediately preceding and following the $i^{\rm
th}$ subcell. In H4, the calculation encompasses both the layers immediately
preceding and following, as well as the virtual layer generated by summing the
energies from two layers downstream and two layers upstream of the $i^{\rm
th}$ subcell. In the S2 arrangement, only the virtual layer computed from the
summation of energies in the directly preceding and following layers is used.

We proceed to assign the updated signal within the current subcell as follows:

$E_{i}=E_{\rm tile}W_{i}/\sum_{j}W_{j}.$ (4.4)

Here, the summation of weights across $j$ serves to normalize the weights,
ensuring that the total signal within the tile remains unchanged, i.e.,
$\sum_{i}E_{i}=E_{\rm tile}$.

The shower’s position is subsequently calculated in a manner analogous to Eq.
4.1:

$\vec{x}_{\rm recon}=\frac{\sum\limits_{i\in\rm
subcells}\vec{x}_{i}w_{i}}{\sum\limits_{i\in\rm subcells}w_{i}}$ (4.5)

where $\vec{x}_{i}$ are the 3d positions of the subcells, and the weights,
$w_{i}$, are determined by

$w_{i}=\max\left(0,w_{0}+\ln{\frac{E_{i}}{E_{\rm tot}}}\right)$ (4.6)

where $E_{i}$ is the energy of the subcell, and $E_{\rm tot}$ is the total
energy of the shower.

This technique extends the principles of the SSA method [11], which focused on
the overlap between vertical and horizontal strips. A notable distinction
between HEXPLIT and the SSA lies in their applicability. The SSA is tailored
for strip geometries in which the layers directly preceding and following a
given layer share identical layouts, e.g. those used in Ref.[9, 10, 22, 23].
On the contrary, HEXPLIT extends this concept to encompass non-strip
geometries and multifold tessellating cycles. In this framework, the
intersection of multiple adjacent layers, each characterized by a distinct
pattern, plays a role in establishing the subcell boundaries and their
associated energies.

For detectors where the cell size varies per layer (as is being considered for
the EIC calorimeter insert design [18]), the weights in Eq. 4.6 should be
adjusted by a factor that is inversely proportional to the area of an
individual cell.

###  4.3 Example reconstruction of a hadronic shower

To illustrate the improvement in positional resolution achieved with HEXPLIT,
we present a section of a simulated neutron shower measured using the H3
configuration in Fig. 2. The top panel displays the energy distribution across
individual cells. In the middle panel, we demonstrate the shower-position
reconstruction using the conventional log-weighting method. In this
representation, only cells that contribute to position determination are shown
with non-zero values. The color spectrum corresponds to the weights used in
determining the shower position, with a threshold of $w_{0}=4.0$ applied. The
bottom panel shows the results when using HEXPLIT. The colors assigned to the
triangular sub-cells represent the weights attributed to the sub-cell
contributions in determining the shower position.

Figure 2:  Top: A segment of a simulated neutron shower example employing the
H3 staggered configuration. Middle: Weights of the hits in the baseline
shower-position reconstruction algorithm. Bottom: Weights of the sub-cells in
the shower-position reconstruction using the HEXPLIT algorithm, incorporating
sub-cell reweighting. The red line depicts the trajectory of the neutron,
extrapolated from its generator-level position and angle. In this coordinate
system, $z=0$ corresponds to the front face of the detector.

##  5 Results

Figure 3 shows the distributions of the radial-position residuals, defined by

$\Delta r=r_{\rm recon}-r_{\rm truth},$ (5.1)

for 50 GeV neutrons in each of hexagonal-cell and square-cell layouts
described in Sec. 2. Here, $r_{\text{recon}}$ is the reconstructed radial
position of the shower, while $r_{\text{truth}}$ represents the true position
of the incoming particle. More specifically, the value of $r_{\rm recon}$ is
given by $\sqrt{x_{\rm recon}^{2}+y_{\rm recon}^{2}}$, where $x_{\rm recon}$
and $y_{\rm recon}$ are given by Eq. 4.1 for the baseline algorithm and Eq.
4.5 for HEXPLIT.

The value of $r_{\text{truth}}$ is analogously defined for the point along the
truth-particle trajectory that aligns with the same longitudinal position as
the reconstructed position (i.e., the log-weighted center of the shower, as
represented by the $z$ component in Eq. 4.1 or Eq. 4.5). This formulation of
position residuals mirrors the approach used in Ref. [20].

For the $w_{0}$ parameter, we used the following functional form

$w_{0}(E)=a+b\,\log\frac{E}{\rm 50\,GeV}+c\,(\log\frac{E}{\rm 50\,GeV})^{2},$
(5.2)

Where $E$ represents the reconstructed energy in GeV, and $a$ and $b$ are
constants. The values of $a$, $b$, and $c$ used in this study are presented in
Table 2. These values were determined through inspection to yield optimal
outcomes.

configuration | $a$ | $b$ | $c$
---|---|---|---
unstaggered hexagons | 5.9 | 0.33 | 0.0
H3 baseline | 4.0 | 0.019 | 0.59
H4 baseline | 4.2 | 0.045 | 0.32
H3 HEXPLIT | 5.0 | 0.89 | 0.17
H4 HEXPLIT | 5.0 | 0.65 | 0.31
Table 2: Values of $a$, $b$, and $c$ used to calculate the $w_{0}$ value used
in reconstruction.

In the top row of Fig. 3, we present the residual distributions obtained for
50 GeV neutrons for five distinct scenarios: unstaggered hexagons, H3
staggering with the baseline reconstruction algorithm, H4 staggering with the
baseline reconstruction algorithm, H3 staggering with HEXPLIT, and H4
staggering with HEXPLIT.

To derive the resolutions, we did fit each distribution with a Gaussian
function and extracted the $\sigma$ values from these fits. Without the
implementation of staggering, the resolution measured at 50 GeV is
10.7$\pm$0.3 mm. By employing H3 staggering, we achieved an improved
resolution of 8.3$\pm$0.4 mm using the baseline reconstruction, which further
improved to 6.0$\pm$0.2 mm with HEXPLIT.

For configurations involving H4 staggering, the results were somewhat better
than those of H3 and unstaggered in the context of baseline reconstruction,
resulting in a resolution of 7.6$\pm$0.5 mm. With the integration of HEXPLIT,
the resolution advanced significantly to 5.1$\pm$0.3 mm.

For completeness, the results for 50 GeV neutrons with square
tessellations—both staggered and unstaggered—are shown in the bottom row of
Fig. 3. Unstaggered squares resulted in a resolution of 10.9$\pm$0.2 mm. For
staggered squares using baseline reconstruction, the resolution was slightly
better but similar, at 8.4$\pm$0.5 mm. On the other hand, staggered squares
employing HEXPLIT reconstruction achieved a resolution of 6.2$\pm$0.3 mm,
representing a 43% improvement over unstaggered squares. The optimal
resolutions were attained by using specific $w_{0}$ values for these square-
cell configurations: 6.3 for unstaggered squares, 4.5 for staggered squares
with baseline reconstruction, and 5.0 for staggered squares with HEXPLIT
reconstruction.

Figure 3: Residual distributions of neutron-shower position reconstruction.
Top row: results with the unstaggered hexagons configuration (left), staggered
configuration H3 (middle) and the staggered configuration H4 (right). Bottom
row: results with the unstaggered squares configuration (left), and the
staggered configuration S2 (right). The results for the staggered
configurations are shown separately with the baseline reconstruction (blue,
open symbols) and with the HEXPLIT algorithm (orange, filled symbols). Each of
these distributions are fit to Gaussian functions (curves).

In the top panel of Fig. 4, we present the position resolution as a function
of neutron energy and compare the results between the unstaggered geometry and
the H3 and H4 staggered configurations. For the staggered layouts, a
comparison is made between results obtained with the baseline reconstruction
and those achieved with HEXPLIT. The bottom panel depicts the extent of
improvement achieved by the staggered layouts when compared to the unstaggered
arrangement.

At the lowest energies, the unstaggered and both staggered layouts have about
the same resolutions when employing the baseline algorithm. With the
application of the HEXPLIT algorithm, there is about a 30% improvement for
both staggered layouts. As the neutron energies increase, the H4 staggering
outperforms the H3 staggering, which in turn exhibits better resolution than
the unstaggered hexagonal configuration.

For mid-to-high energies, the H4 staggering with the baseline reconstruction
results in a substantial improvement, up to 50% (depending on the energy) over
the unstaggered configuration. With the HEXPLIT algorithm, the relative
improvement of H4 in comparison to the unstaggered configuration is even more
pronounced, from about 35% to 60% across the full energy range.

Figure 4: Top Panel: Neutron-shower position resolution comparisons among
different configurations, including the unstaggered hexagon configuration
(black), H3 staggered configuration (blue), and the H4 staggered configuration
(orange). For the staggered simulations, both resolutions obtained with (solid
symbols) and without (open symbols) HEXPLIT are displayed. The horizontal line
represents the resolution anticipated from the side length divided by
$\sqrt{12}$ for squares with the same area as the hexagons used in the
simulations. Bottom Panel: Ratio depicting the resolutions in the staggered
hexagonal configurations relative to those of the unstaggered hexagon
configuration.

To investigate the feasibility of the staggered layouts and the HEXPLIT
algorithm in electromagnetic calorimetry, we repeated our studies for single-
photon events, as described in more detail in Appendix A. We found that the
position resolutions obtained for photon events with the baseline
reconstruction improved substantially with both the H3 and H4 staggering
layouts. However, the HEXPLIT algorithm did not further enhance the
performance as it did for neutron showers. Overall, the resolutions were
improved twofold through the use of H4, similar to the neutron results.

We explored the impact of cell size on resolution by conducting 50 GeV neutron
simulations with hexagons that were either twice as large or half as large as
the original cells. The resolutions for all three sizes are presented in Fig.
5. Our findings indicate that in every case, staggered designs outperform the
unstaggered design. Additionally, H4 staggering produces better outcomes than
H3, and both configurations show improvement when HEXPLIT is applied. We also
observed that resolution improvement reaches a point of diminishing returns as
the cell sizes decrease.

Figure 5: Neutron-shower position resolution comparisons among different
configurations and cell sizes, including the unstaggered hexagon configuration
(black), H3 staggered configuration (blue), and the H4 staggered configuration
(orange). For the staggered simulations, both resolutions obtained with (solid
symbols) and without (open symbols) HEXPLIT are displayed. Bottom Panel: Ratio
depicting the resolutions in the staggered hexagonal configurations relative
to those of the unstaggered hexagon configuration.

##  6 Summary and Conclusions

We have developed an algorithm and introduced innovative tessellation patterns
that significantly enhance position resolution for hadronic showers detected
by high-granularity calorimeters. Our findings underscore the performance
improvements achieved by incorporating staggered hexagonal cells in cyclic
tessellating patterns, resulting in enhancements of up to nearly double the
performance compared to non-staggered configurations.

These findings highlight the substantial potential of staggered
configurations, which can also be harnessed in machine-learning applications
[24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38], particularly
those capable of accommodating complex geometries [27, 39, 40, 41, 42, 43].

The improved position resolution has the potential to enhance the performance
of the particle-flow algorithm [2] for a given detector granularity.
Alternatively, this technique can be employed to optimize costs by reducing
the necessary number of channels in the design of high-granularity
calorimeters, all while maintaining the desired level of spatial resolution.

These outcomes pave the way for future explorations into the viability of
employing staggered scintillator tiling in other ongoing high-granularity
sampling calorimeter projects, including those at the HL-LHC [3], ILC [44],
CEPC [45], and the EIC [18, 46].

##  7 Code Availability

An implementation of the HEXPLIT algorithm, along with example simulated
files, can be found in Ref. [47].

## Acknowledgments

We thank Sean Preins for insightful discussions on tessellation with hexagons
in the design of the calorimeter insert. We thank members of the California
EIC consortium for their feedback on our work. This work was supported by MRPI
program of the University of California Office of the President, award number
00010100. S.P also acknowledges support from the Jefferson Lab EIC Center
Fellowship. M.A acknowledges support through DOE Contract No. DE-
AC05-06OR23177 under which Jefferson Science Associates, LLC operates the
Thomas Jefferson National Accelerator Facility.

## References

  * [1] F. Sefkow, A. White, K. Kawagoe, R. Pöschl, and J. Repond, “Experimental tests of particle flow calorimetry,” Rev. Mod. Phys. 88 (Feb, 2016) 015003. https://link.aps.org/doi/10.1103/RevModPhys.88.015003.
  * [2] M. Thomson, “Particle flow calorimetry and the PandoraPFA algorithm,” Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 611 no. 1, (Nov, 2009) 25–40. https://doi.org/10.1016%2Fj.nima.2009.09.009.
  * [3] CMS Collaboration, “The Phase-2 Upgrade of the CMS Endcap Calorimeter,” tech. rep., CERN, Geneva, 2017.  https://cds.cern.ch/record/2293646.
  * [4] M. N. on behalf of the CMS Collaboration, “The CMS HGCAL silicon region architecture specification and optimisation,” JINST 17 no. 03, (Mar, 2022) C03010. https://dx.doi.org/10.1088/1748-0221/17/03/C03010.
  * [5] ATLAS Collaboration, “ATLAS tile calorimeter: Technical design report,”.
  * [6] CALICE Collaboration, J. Repond et al., “Design and electronics commissioning of the physics prototype of a Si-W electromagnetic calorimeter for the International Linear Collider,” JINST 3 no. 08, (Aug, 2008) P08001. https://dx.doi.org/10.1088/1748-0221/3/08/P08001.
  * [7] ILD Concept Group Collaboration, H. Abramowicz et al., “International Large Detector: Interim Design Report,” arXiv:2003.01116 [physics.ins-det].
  * [8] CALICE Collaboration, “Study of position and angular resolution for electron showers measured with the electromagnetic SiW prototype,” tech. rep., 2009.  http://cds.cern.ch/record/2669469.
  * [9] CALICE Collaboration, K. Francis et al., “Performance of the first prototype of the CALICE scintillator strip electromagnetic calorimeter,” Nucl. Instrum. Meth. A 763 (2014) 278–289, arXiv:1311.3761 [physics.ins-det].
  * [10] CALICE Collaboration, J. Repond et al., “Construction and Response of a Highly Granular Scintillator-based Electromagnetic Calorimeter,” Nucl. Instrum. Meth. A 887 (2018) 150–168, arXiv:1707.07126 [physics.ins-det].
  * [11] K. Kotera, D. Jeans, A. Miyamoto, and T. Takeshita, “A novel strip energy splitting algorithm for the fine granular readout of a scintillator strip electromagnetic calorimeter,” Nucl. Instrum. and Meth. A 789 (2015) 158–164. https://www.sciencedirect.com/science/article/pii/S0168900215004544.
  * [12] L. M. S. de Silva and F. Simon, “Effects of misalignment on response uniformity of SiPM-on-tile technology for highly granular calorimeters,” JINST 15 no. 06, (2020) P06030, arXiv:2004.05066 [physics.ins-det].
  * [13] A. Belloni et al., “Test beam study of SiPM-on-tile configurations,” JINST 16 no. 07, (2021) P07022, arXiv:2102.08499 [physics.ins-det].
  * [14] CALICE Collaboration, “Design, Construction and Commissioning of a Technological Prototype of a Highly Granular SiPM-on-tile Scintillator-Steel Hadronic Calorimeter,” arXiv:2209.15327 [physics.ins-det].
  * [15] M. Arratia, L. Garabito Ruiz, J. Huang, S. J. Paul, S. Preins, and M. Rodriguez, “Studies of time resolution, light yield, and crosstalk using SiPM-on-tile calorimetry for the future Electron-Ion Collider,” JINST 18 no. 05, (2023) P05045, arXiv:2302.03646 [physics.ins-det].
  * [16] M. Frank, F. Gaede, C. Grefe, and P. Mato, “DD4hep: A Detector Description Toolkit for High Energy Physics Experiments,” J. Phys. Conf. Ser. 513 (2014) 022010.
  * [17] GEANT4 Collaboration, S. Agostinelli et al., “GEANT4–a simulation toolkit,” Nucl. Instrum. Meth. A 506 (2003) 250–303.
  * [18] M. Arratia et al., “A high-granularity calorimeter insert based on SiPM-on-tile technology at the future Electron-Ion Collider,” Nucl. Instrum. Meth. A 1047 (2023) 167866, arXiv:2208.05472 [physics.ins-det].
  * [19] Sebouh J. Paul, “Staggered Segmentations plugin for DD4HEP.”  https://github.com/sebouh137/staggered_segmentations.
  * [20] CMS HGCAL Collaboration, B. Acar et al., “Response of a CMS HGCAL silicon-pad electromagnetic calorimeter prototype to 20–300 GeV positrons,” JINST 17 no. 05, (May, 2022) P05022. https://dx.doi.org/10.1088/1748-0221/17/05/P05022.
  * [21] N. Akchurin et al., “First beam tests of prototype silicon modules for the CMS High Granularity Endcap Calorimeter,” JINST 13 no. 10, (Oct, 2018) P10023. https://dx.doi.org/10.1088/1748-0221/13/10/P10023.
  * [22] N. Tsuji et al., “Highly granular scintillator-strip electromagnetic calorimeter for future higgs factories,” Nucl. Instrum. Meth. A 1047 (2023) 167861. https://www.sciencedirect.com/science/article/pii/S0168900222011536.
  * [23] M. Arratia, R. Milton, S. J. Paul, B. Schmookler, and W. Zhang, “A Few-Degree Calorimeter for the future Electron-Ion Collider,” arXiv:2307.12531 [physics.ins-det].
  * [24] M. Paganini, L. de Oliveira, and B. Nachman, “CaloGAN : Simulating 3D high energy particle showers in multilayer electromagnetic calorimeters with generative adversarial networks,” Phys. Rev. D 97 no. 1, (2018) 014021, arXiv:1712.10321 [hep-ex].
  * [25] M. Paganini, L. de Oliveira, and B. Nachman, “Accelerating Science with Generative Adversarial Networks: An Application to 3D Particle Showers in Multilayer Calorimeters,” Phys. Rev. Lett. 120 no. 4, (2018) 042003, arXiv:1705.02355 [hep-ex].
  * [26] D. Belayneh et al., “Calorimetry with deep learning: particle simulation and reconstruction for collider physics,” Eur. Phys. J. C 80 no. 7, (2020) 688, arXiv:1912.06794 [physics.ins-det].
  * [27] S. R. Qasim, J. Kieseler, Y. Iiyama, and M. Pierini, “Learning representations of irregular particle-detector geometry with distance-weighted graph networks,” Eur. Phys. J. C 79 no. 7, (2019) 608, arXiv:1902.07987 [physics.data-an].
  * [28] F. A. Di Bello, S. Ganguly, E. Gross, M. Kado, M. Pitt, L. Santi, and J. Shlomi, “Towards a Computer Vision Particle Flow,” Eur. Phys. J. C 81 no. 2, (2021) 107, arXiv:2003.08863 [physics.data-an].
  * [29] E. Buhmann, S. Diefenbacher, E. Eren, F. Gaede, G. Kasieczka, A. Korol, and K. Krüger, “Getting High: High Fidelity Simulation of High Granularity Calorimeters with High Speed,” Comput. Softw. Big Sci. 5 no. 1, (2021) 13, arXiv:2005.05334 [physics.ins-det].
  * [30] N. Akchurin, C. Cowden, J. Damgov, A. Hussain, and S. Kunori, “On the use of neural networks for energy reconstruction in high-granularity calorimeters,” JINST 16 no. 12, (2021) P12036, arXiv:2107.10207 [physics.ins-det].
  * [31] J. Pata, J. Duarte, J.-R. Vlimant, M. Pierini, and M. Spiropulu, “MLPF: Efficient machine-learned particle-flow reconstruction using graph neural networks,” Eur. Phys. J. C 81 no. 5, (2021) 381, arXiv:2101.08578 [physics.data-an].
  * [32] C. Neubüser, J. Kieseler, and P. Lujan, “Optimising longitudinal and lateral calorimeter granularity for software compensation in hadronic showers using deep neural networks,” Eur. Phys. J. C 82 no. 1, (2022) 92, arXiv:2101.08150 [physics.ins-det].
  * [33] N. Akchurin, C. Cowden, J. Damgov, A. Hussain, and S. Kunori, “Perspectives on the Calibration of CNN Energy Reconstruction in Highly Granular Calorimeters,” arXiv:2108.10963 [physics.ins-det].
  * [34] E. Buhmann, S. Diefenbacher, D. Hundhausen, G. Kasieczka, W. Korcari, E. Eren, F. Gaede, K. Krüger, P. McKeown, and L. Rustige, “Hadrons, better, faster, stronger,” Mach. Learn. Sci. Tech. 3 no. 2, (2022) 025014, arXiv:2112.09709 [physics.ins-det].
  * [35] G. R. Khattak, S. Vallecorsa, F. Carminati, and G. M. Khan, “Fast simulation of a high granularity calorimeter by generative adversarial networks,” Eur. Phys. J. C 82 no. 4, (2022) 386, arXiv:2109.07388 [physics.ins-det].
  * [36] M. Chadeeva and S. Korpachev, “Machine-learning-based prediction of parameters of secondaries in hadronic showers using calorimetric observables,” JINST 17 no. 10, (2022) P10031, arXiv:2205.12534 [physics.ins-det].
  * [37] S. R. Qasim, N. Chernyavskaya, J. Kieseler, K. Long, O. Viazlo, M. Pierini, and R. Nawaz, “End-to-end multi-particle reconstruction in high occupancy imaging calorimeters with graph neural networks,” Eur. Phys. J. C 82 no. 8, (2022) 753, arXiv:2204.01681 [physics.ins-det].
  * [38] V. Mikuni and B. Nachman, “Score-based generative models for calorimeter shower simulation,” Phys. Rev. D 106 no. 9, (2022) 092009, arXiv:2206.11898 [hep-ph].
  * [39] E. Buhmann, S. Diefenbacher, E. Eren, F. Gaede, G. Kasieczka, A. Korol, W. Korcari, K. Krüger, and P. McKeown, “CaloClouds: Fast Geometry-Independent Highly-Granular Calorimeter Simulation,” arXiv:2305.04847 [physics.ins-det].
  * [40] J. Liu, A. Ghosh, D. Smith, P. Baldi, and D. Whiteson, “Generalizing to new geometries with Geometry-Aware Autoregressive Models (GAAMs) for fast calorimeter simulation,” arXiv:2305.11531 [physics.ins-det].
  * [41] F. T. Acosta, V. Mikuni, B. Nachman, M. Arratia, K. Barish, B. Karki, R. Milton, P. Karande, and A. Angerami, “Comparison of Point Cloud and Image-based Models for Calorimeter Fast Simulation,” arXiv:2307.04780 [cs.LG].
  * [42] O. Amram and K. Pedro, “CaloDiffusion with GLaM for High Fidelity Calorimeter Simulation,” arXiv:2308.03876 [physics.ins-det].
  * [43] H. Hashemi, N. Hartmann, S. Sharifzadeh, J. Kahn, and T. Kuhr, “Ultra-High-Resolution Detector Simulation with Intra-Event Aware GAN and Self-Supervised Relational Reasoning,” arXiv:2303.08046 [physics.ins-det].
  * [44] T. Behnke, J. E. Brau, P. N. Burrows, J. Fuster, M. Peskin, M. Stanitzki, Y. Sugimoto, S. Yamada, and H. Yamamoto, “The international linear collider technical design report - volume 4: Detectors,” 2013.
  * [45] CEPC Study Group Collaboration, M. Dong et al., “CEPC Conceptual Design Report: Volume 2 - Physics & Detector,” arXiv:1811.10545 [hep-ex].
  * [46] F. Bock et al., “Design and Simulated Performance of Calorimetry Systems for the ECCE Detector at the Electron Ion Collider,” Nucl. Instrum. Meth. A 1055 (2023) 168464, arXiv:2207.09437 [physics.ins-det].
  * [47] Sebouh J. Paul, “Example files and Python code for the HEXPLIT algorithm,” 2023.  https://doi.org/10.5281/zenodo.8247336.

##  Appendix A Electromagnetic calorimetry

We repeated our studies with photons instead of neutrons. We used a separate
optimization of the $w_{0}$ parameters for this study, using the parameters in
Table 3. We show the distributions of the position residuals for 50 GeV
photons in Fig. 6, analogous to Fig. 3. We found that for each of the layouts
considered (unstaggered hexagons, H3, and H4), the position resolutions
obtained using the baseline reconstruction for photons are several times
better than those for neutrons (see Fig. 3). Furthermore, the staggering
pattern H3 offers better resolution than the unstaggered hexagons, and the H4
staggering has better resolution than H3. However, for both the H3 and H4
staggering patterns, HEXPLIT is found to yield performance similar to the
baseline algorithm. This may be a limitation of the detector geometry and the
cell size considered.

We show the dependence of the position resolution on the photon energy in Fig.
7, analogous to Fig. 4. The resolution with the H3 staggering using the
baseline algorithm is around 30% better than that of the unstaggered hexagons
across the full energy range, while the H4 staggering is about 50% better than
the unstaggered hexagons.

configuration | $a$ | $b$ | $c$
---|---|---|---
unstaggered hexagons | 7.3 |  $-$0.14 |  $-$0.088
H3 baseline | 6.3 | 0.25 |  $-$0.23
H4 baseline | 5.4 | 0.37 | 0.0
H3 HEXPLIT | 7.6 | 0.38 | 0.0
H4 HEXPLIT | 6.0 | -0.46 | 0.099
Table 3: Values of $a$, $b$, and $c$ used to calculate the $w_{0}$ value used
in the reconstruction of photons. Figure 6: Residual distributions of photon-
shower position reconstruction. Results with the unstaggered hexagons
configuration (left), staggered configuration H3 (middle) and the staggered
configuration H4 (right). The results for the staggered configurations are
shown separately with the baseline reconstruction (blue, open symbols) and
with the HEXPLIT algorithm (orange, filled symbols). Each of these
distributions are fit to Gaussian functions (curves).  Figure 7: Photon-shower
position resolution comparisons among different configurations, including the
unstaggered hexagon configuration (black), H3 staggered configuration (blue),
and the H4 staggered configuration (orange). For the staggered simulations,
both resolutions obtained with (solid symbols) and without (open symbols)
HEXPLIT are displayed. Bottom Panel: Ratio depicting the resolutions in the
staggered hexagonal configurations relative to those of the unstaggered
hexagon configuration.