#  Localization and slow-thermalization in a cluster spin model

Yoshihito Kuno1 Takahiro Orito2 Ikuo Ichinose3 1Department of Physics,
Graduate School of Science, Tsukuba University, Tsukuba, Ibaraki 305-8571,
Japan  2Graduate School of Advanced Science and Engineering, Hiroshima
University, 739-8530, Japan  3Department of Applied Physics, Nagoya Institute
of Technology, Nagoya, 466-8555, Japan

###### Abstract

New cluster spin model with interactions and disorder is studied. In certain
type of interactions, we find an extensive number of local integrals of motion
(LIOMs), which are a modified version of the stabilizers, but not dressed
spins. These LIOMs can be defined for any strength of the interactions and
disorder. Hence, even under the presence of interactions, integrability is
held, and all eigenvalues are labeled by these LIOMs. Localization is, then,
expected to occur. We numerically investigate this expectation in detail. In
particular, ergodicity breaking dynamics is demonstrated where symmetry-
protected-topological order also preserves for specific initial states.

Introduction.— The notions of localization and thermalization are deeply
related Nandkishore2015 ; Abanin . In a closed system separated from
environment, if states are localized, the system does not thermalize. This
character persists even in interacting many-body systems, phenomenon called
many-body-localization (MBL) Basko2006 ; Bardarson . Such non- or slow-
thermalization dynamics is a novel example for the breaking of the eigenstate
thermalization hypothesis Rigol2007 ; Nandkishore2015 . Also, weak
localization and weak ergodicity breaking in some specific models such as
quantum scar states are now getting a lot of attention Bernien ; Turner2018 ;
Choi2019 ; Ho2019 ; Iadecola2020 ; Bluvstein2020 ; McClarty2020 ; Shibata2020
; Lee2020 ; KMH2020 ; Serbyn2021 ; Jeyaretnam ; KMH2020 ; KMH2021 . What types
of models exhibit localization or MBL phenomena, i.e., non- or slow-
thermalization dynamics, and what constraints induce them are significant
questions. To answer them, a key concept is an extensive number of local
integrals of motions (LIOMs) Serbyn2013 ; Huse2014 ; Imbrie2016 ; Imbrie2017 .
Investigating the LIOMs in specific models is important and useful to
scrutinize localization phenomena from the general perspective. Constructing
some classification schema of localization from a bottom-up approach by the
study on various concrete examples is an important subject Serbyn2021 .

In this Letter, we study a specific example of systems, in which LIOMs are
obtained explicitly. The LIOMs are compact support in the system with certain
interactions and disorder, and the system exhibits localization phenomena in
its dynamics. Model describing this system is an extended version of the
cluster model Briegel2001 ; Pachos2004 ; Son ; Smacchia , which includes
interactions and disorder. In the pure cluster model, the model exhibits
localization phenomena with symmetry-protected-topological (SPT) order
Bauer2013 ; Bahri ; Vasseur2016 ; Parameswaran ; Decker ; Kuno2019 ;
Wahl2020_1 ; Wahl2020_2 ; Wahl2020_3 ; Kemp2020 ; Sahay2021 ; Duque2021 . For
the extended version of the cluster model, we find localization phenomena such
as slow-thermalization analog to Anderson localization. The extensive number
of LIOMs are modified counterparts of the stabilizers in the pure cluster
model, and they label all eigenstates of the interacting and disordered
Hamiltonian. The compact support of the LIOMs is related to the fact that the
interactions only locally mix many-body states. Expected localization
phenomena in the model including ergodicity breaking are numerically
demonstrated. We, furthermore, investigate how localization and the SPT order
are destroyed as an additional interaction is added, although they are quite
stable regardless of the strength of the original interactions and disorder.

Model.— In this Letter, we work on an extended version of the cluster spin
model with interactions and disorder, the original model of which was studied
in Ref. Briegel2001 ; Pachos2004 ; Son ; Smacchia ,

$\displaystyle
H=\sum^{L/2-1}_{\ell=0}[J_{2\ell}K_{2\ell}+J_{2\ell+1}K_{2\ell+1}]+H_{\rm
int}$ (1) $\displaystyle H_{\rm int}=\sum^{L/2-1}_{\ell=0}gV_{\ell},$ (2)

where $K_{j}=\sigma^{z}_{j-1}\sigma^{x}_{j}\sigma^{z}_{j+1}$ is a stabilizer
operator composed of Pauli matrices,
$V_{\ell}=N_{\ell-1}[K^{+}_{2\ell}K^{-}_{2\ell+1}+K^{-}_{2\ell}K^{+}_{2\ell+1}]$,
$N_{\ell}=K_{2\ell}+K_{2\ell+1}$ and
$K^{\pm}_{j}=\frac{1}{2}(\sigma^{z}_{j}\mp
i\sigma^{z}_{j-1}\sigma^{y}_{j}\sigma^{z}_{j+1})$. $K^{\pm}_{j}$ is a raising
and lowering operator for cluster state Jeyaretnam , $J_{j}$ is site-dependent
potential of the stabilizer $K_{j}$ and $g$ is strength of the interactions.
The stabilizer $K_{j}$ is a dressed spin, which satisfies $K^{2}_{j}=1$ and
SU(2)-spin algebra, $[K_{j},K^{\pm}_{j}]=\pm 2K_{j}$. The Hamiltonian $H$ has
$\mathbb{Z}_{2}\times\mathbb{Z}^{T}_{2}$ symmetry, which is composed of the
global spin flip $\prod^{L-1}_{j=0}\sigma^{x}_{j}$ and the complex conjugation
(time-reversal operation) Verresen2017 ; Verresen2018 ; Smith2020 . For $g=0$,
the set of the stabilizer $\\{K_{j}\\}$ becomes an extensive number of LIOMs
since $[H_{g=0},K_{j}]=0$ and $[K_{j},K_{j^{\prime}}]=0$ for any $j$ and
$j^{\prime}$. When the site-dependent coupling $J_{j}$ ($j=0,1,\cdots,L-1$) is
randomly valued, the system exhibits localization with short-range
entanglement coming from the cluster state, where every eigenstate is uniquely
labeled by the set of eigenvalues of the LIOMs, $K_{j}=\pm 1$. This
localization is stable againt perturbations (interactions) preserving the
symmetry, namely, topological MBL, which has been studied extensively Bahri ;
Decker ; Vasseur2016 ; Parameswaran ; Decker ; Kuno2019 ; Wahl2020_1 ;
Wahl2020_2 ; Wahl2020_3 ; Kemp2020 ; Sahay2021 ; Duque2021 . In this Letter,
we focus on $g\neq 0$ and a type of disorder such as
$J_{2\ell}=-J_{2\ell+1}=\lambda_{\ell}$, where $\lambda_{\ell}$’s are
uniformly-distributed random variables, $\lambda_{\ell}\in[-W,W]$. In this
case, the system has modified LIOMs as we show shortly.

Figure 1: Schematic figure of the Model, Eq. (3). The upper chain is the
original spin-lattice. The original spin-lattice is mapped onto the cluster-
basis lattice in the lower chain. The green vertical dashed line represents
the entanglement cut for a pair of subsystems in the calculation of the EE.
Figure 2: Identification of eigenstates and their LIOM-eigenvalues for three
unit-cell system ($L=6$). (a) $|\Delta\psi_{k}|^{2}$ for each eigenstate for
$g=0$ case. (b) Each eigenvalue of $K^{a}_{\ell_{0}}$ for $g=0$ case. (c)
$|\Delta\psi_{k}|^{2}$ for each eigenstate for finite $g$ in terms of
$K^{a}_{\ell_{0}}$. (d) Each eigenvalue
$\langle\psi_{k}|K^{a}_{\ell_{0}}|\psi_{k}\rangle$ for finite $g$ case. (e)
$|\Delta\psi_{k}|^{2}$ for each eigenstate for finite $g$ in terms of
$\tilde{K}^{a}_{\ell_{0}}$. (f) Each eigenvalue of $\tilde{K}^{a}_{\ell_{0}}$
for finite $g$ case. All results are obtained for a single-shot disorder.
Energy eigenstates are numbered in the ascendant order.

Modified local integrals of motion.— Under the disorder
$\\{\lambda_{\ell}\\}$, the target cluster model is rewritten as

$\displaystyle H=H_{0}+H_{\rm
int},\>\>H_{0}=\sum^{L/2-1}_{\ell=0}\lambda_{\ell}[K^{a}_{\ell}-K^{b}_{\ell}],$
(3)

where we have introduced unit-cell including two sites $2\ell$ and $2\ell+1$,
and “$a$” and “$b$” indices for even and odd sites in the unit-cell. The
operators $K_{2\ell(2\ell+1)}$ are relabeled as $K^{a(b)}_{\ell}$. The
schematic lattice structure is shown in Fig. 1. Similar structure to
$K^{a(b)}_{\ell}$ was used in other lattice models, e.g., compact localized
states KOI2020 ; OKI2020 ; OKI2021 .

By following our previous study on the Creutz ladder OKI2020 , we find that
there are operators commutative with the total Hamiltonian $H$, which are
given by

$\displaystyle\tilde{K}^{a}_{\ell}=K^{a}_{\ell}+\frac{g}{2\lambda_{\ell}}V_{\ell},\>\>\tilde{K}^{b}_{\ell}=K^{b}_{\ell}-\frac{g}{2\lambda_{\ell}}V_{\ell},$
(4)

and
$H=\sum^{L/2-1}_{\ell=0}\lambda_{\ell}[\tilde{K}^{a}_{\ell}-\tilde{K}^{b}_{\ell}]$.
These operators are compact-support form and satisfy
$[\tilde{K}^{a}_{\ell},\tilde{K}^{a(b)}_{\ell^{\prime}}]=0$ for any $\ell$ and
$\ell^{\prime}$, $[H,\tilde{K}^{a(b)}_{\ell}]=0$. That is, these operators are
LIOMs, which are a modified version of the original stabilizer $K_{j}$. Also,
note that the LIOMs, $\tilde{K}^{a(b)}_{\ell}$, can be defined for any
disorder realizations and interaction strength, $\\{\lambda_{\ell}\\}$ and
$g$, whereas $(\tilde{K}^{a(b)}_{\ell})^{2}\neq 1$, that is, the eigenvalues
of $\tilde{K}^{a(b)}_{\ell}$ are not a discrete value $\pm 1$, but fractional
in general, depending on the values of $g$ and $\lambda_{\ell}$. Hence,
$\tilde{K}^{a(b)}_{\ell}|\psi_{k}\rangle=I^{k}_{a(b),\ell}|\psi_{k}\rangle$,
where $I^{k}_{a(b),\ell}\in\mathbb{R}$ and $|\psi_{k}\rangle$ is $k$-th
eigenstate of $H$. [We use the ascendant order for the energy eigenstates.]
Here, we remark that we can extend systems to the one with long-range
interactions with modified LIOMs beyond compact support, which we comment on
in conclusion.

In the following, we show how the LIOMs label all eigenstates. Since extensive
number of the LIOMs [${\tilde{K}}^{a(b)}_{\ell}$] exist, the system can
exhibit localization nature, somewhat analog to Anderson localization of
cluster-spin basis. Furthermore, as characteristic dynamical properties, non-
or slow-thermalized dynamics appears Rigol2007 .

Study of small system.— We study a small system with $L=6$ to see how the
system is affected by the interactions, how the eigenstates are characterized
by LIOMs and how localization properties emerge.

In this small system, we first focus on the $\ell$-th unit-cell in the
stabilizer lattice in Fig. 1 and investigate the eigenvalues and eigenstates
of the non-interacting Hamiltonian $H_{0}$. For any $\lambda_{\ell}$, there
are three energy levels, $\epsilon=-2\lambda_{\ell}$, $0$, $2\lambda_{\ell}$,
where their eigenstates are given by $|\bar{1}1\rangle_{\ell}$ for
$\epsilon=-2\lambda_{j}$, two orthogonal linear-superposed states denoted by
$\alpha|11\rangle_{\ell}+\beta|\bar{1}\bar{1}\rangle_{\ell}$ for $\epsilon=0$,
and $|1\bar{1}\rangle$ for $\epsilon=2\lambda_{j}$, where we have introduced a
notation such as
$|1\bar{1}\rangle_{\ell}=K^{a+}_{\ell}K^{b-}_{\ell}|\Uparrow\rangle_{\ell}$
($|\Uparrow\rangle_{\ell}$ is all up states around $\ell$ unit-cells in the
original spin lattice). Therefore,
$K^{a}_{\ell}|1{\bar{1}}\rangle_{\ell}=|1\bar{1}\rangle_{\ell}$,
$K^{b}_{\ell}|1{\bar{1}}\rangle_{\ell}=-|1{\bar{1}}\rangle_{\ell}$. The
presence of an arbitrary linear-superposed state at zero energy of $H_{0}$
prevents unique-labeling of eigenstates by the LIOMs. A pair of the doubly
degenerate states can be chosen arbitrarily as long as they are orthogonal to
each other. Thus, we take $|\bar{1}\bar{1}\rangle_{\ell}$ and
$|11\rangle_{\ell}$ as a pair of orthogonal eigenstates with zero energy of
$H_{0}$. To impose this choice in practical calculations, we introduce very
small random potential, $\sum_{\ell}\delta
h_{\ell}[K^{a}_{\ell}+K^{b}_{\ell}]$ with $\delta h_{\ell}\in[-\delta h,\delta
h]$, $\delta h=0.5\times 10^{-5}W$. This ‘fictitious’ disorder gives little
effect to the entire physics, especially, to the dynamical behavior of the
system. Under this manipulation, “cluster-basis” eigenstates of $H_{0}$ are
described as
$|\psi(\\{a_{\ell}\\},\\{b_{\ell}\\})\rangle=\prod^{L/2}_{\ell=0}K^{a\>a_{\ell}}_{\ell}K^{b\>b_{\ell}}_{\ell}|\Uparrow\rangle$
Jeyaretnam , where $L$ is an even number, the sets of
$\\{a_{\ell}\\},\\{b_{\ell}\\}$ are a sequence of $+$ and $-$ labels and
$|\Uparrow\rangle$ is a $L$-site ferromaginetic state with all spins up in the
original spin basis.

Removing the degeneracy of the zero-energy state of $H_{0}$ in single unit-
cell as explained in the above, we study effects of the interaction $H_{\rm
int}$ on the cluster-basis eigenstates of $H_{0}$. Here, an interaction term
$V_{\ell}$ acts over two unit-cells as shown in the lower lattice in Fig. 1.
Only four cluster-basis states on the $\ell-1$ and $\ell$ unit-cells are
changed: (i) $|11\rangle_{\ell-1}|1\bar{1}\rangle_{\ell}$, (ii)
$|11\rangle_{\ell-1}|\bar{1}1\rangle_{\ell}$, (iii)
$|\bar{1}\bar{1}\rangle_{\ell-1}|1\bar{1}\rangle_{\ell}$, (iv)
$|\bar{1}\bar{1}\rangle_{\ell-1}|\bar{1}1\rangle_{\ell}$. The interaction
$V_{\ell}$ mixes (i) and (ii) ((iii) and (iv)), and creates a mixing state of
$|1\bar{1}\rangle_{\ell}$ and $|\bar{1}1\rangle_{\ell}$. For the other
cluster-basis states on the two unit-cells are a null state of $V_{\ell}$.

From the above observation about the action of $V_{\ell}$, characteristic
eigenstates for the $L=6$ interacting system of $H$ are obtained
straightforwardly. As an examples,
$|\psi^{L=6}_{1,\pm}\rangle=|11\rangle_{0}|s^{\pm}\rangle_{1}|11\rangle_{2}$
where
$|s^{\pm}\rangle_{1}=\alpha^{\pm}_{1}|1\bar{1}\rangle_{1}+\alpha^{\pm}_{2}|\bar{1}1\rangle_{1}$
and
$(\alpha^{\pm}_{1},\alpha^{\pm}_{2})=[4g^{2}+(2\lambda_{1}-\epsilon^{\pm}_{1})^{2}]^{-1/2}(2g,\epsilon^{\pm}_{1}-2\lambda_{1})$
with $\epsilon^{\pm}_{1}=\pm 2[\lambda^{2}_{\ell}+g^{2}]^{1/2}$. The state,
$|\psi^{L=6}_{1,\pm}\rangle$, is an eigenstate for all LIOMs,
$\tilde{K}^{a(b)}_{\ell}$, with integer or fractional eigenvalues, e.g.,
$\tilde{K}^{a}_{1}|\psi^{L=6}_{1,\pm}\rangle=\frac{\epsilon^{\pm}_{1}}{2\lambda_{1}}|\psi^{L=6}_{1,\pm}\rangle$
and
$\tilde{K}^{b}_{1}|\psi^{L=6}_{1,\pm}\rangle=-\frac{\epsilon^{\pm}_{1}}{2\lambda_{1}}|\psi^{L=6}_{1,\pm}\rangle$.
From this observation of $|\psi^{L=6}_{1,\pm}\rangle$, certain cluster-basis
eigenstates of $H_{0}$, $|\psi(\\{a_{\ell}\\},\\{b_{\ell}\\})\rangle$ are
mixed by the interactions $H_{\rm int}$, however, the mixing is local and only
small numbers of cluster-basis eigenstates are affected. These observations
indicate that global hybridization does not occur by the interactions $H_{\rm
int}$ due to the presence of the extensive number of the modified LIOMs,
$\tilde{K}^{a(b)}_{\ell}$. [We give observation of the general structure of
the Hilbert space for the $L$-site system in Supp0 .]

Following the above analytical observation, we numerically study the LIOMs in
the $L=6$ system Quspin . [We employ the periodic boundary condition.] We
observe whether the operators $\\{\tilde{K}^{a}_{j}\\}$ operate as the LIOMs
for all eigenstates, i.e., all energy eigenstates are the eigenstates of the
LIOMs. To this end, we first examine whether the resulting states obtained by
acting $\tilde{K}^{a}_{\ell_{0}}$ on energy eigenstates are proposional to the
original ones. Numerically, for $k$-th eigenstate $|\psi_{k}\rangle$, we
define
$|\tilde{\psi}_{k}\rangle\equiv\tilde{K}^{a}_{\ell_{0}}|\psi_{k}\rangle/\sqrt{|\tilde{K}^{a}_{\ell_{0}}|\psi_{k}\rangle|}$
and calculate
$|\Delta\psi_{k}|^{2}=||\psi_{k}\rangle-|\tilde{\psi}_{k}\rangle|^{2}$. Then,
if $|\Delta\psi_{k}|^{2}=0$ or $4$, $|\psi_{k}\rangle$ is also an eigenstate
of $\tilde{K}^{a}_{\ell_{0}}$ with the eigenvalue,
$I^{k}_{a,\ell_{0}}=\langle\psi_{k}|\tilde{K}^{a}_{\ell_{0}}|\psi_{k}\rangle$.
Here, we set $\ell_{0}=1$ for the practical calculation.

We show numerical results. For the $g=0$ case, all eigenstates for single-shot
disorder realization are labeled by the original stabilizers,
$\\{K^{a(b)}_{\ell}\\}$. Each $|\Delta\psi_{k}|^{2}$ takes $0$ or $4$ and the
LIOMs’ eigenvalues are $I^{k}_{a,1}=\pm 1$ as shown in Figs. 2 (a) and (b).
For the $g=1$ case, we first operate the original $\\{K^{a(b)}_{\ell}\\}$ to
energy eigenstates for $g=1$ to see that $\\{K^{a(b)}_{\ell}\\}$ are not
genuine LIOMs. As shown in Figs. 2 (c) and (d), some of $|\Delta\psi_{k}|^{2}$
deviate from $0/4$ and
$\langle\psi_{k}|\tilde{K}^{a}_{\ell_{0}}|\psi_{k}\rangle$ deviate from $\pm
1$. On the other hand for the operation of $\tilde{K}^{a(b)}_{\ell_{0}}$,
$|\Delta\psi_{k}|^{2}$ takes $0$ or $4$, that is, $\tilde{K}^{a}_{\ell_{0}}$
is a genuine LIOM with eigenvalues $\\{I^{k}_{a,1}\\}$, where some of
$I^{k}_{a,1}$ takes a disorder-dependent fractional value. Also, note that the
interactions $H_{\rm int}$ has some large kernel space, $H_{\rm
int}|\psi_{k}\rangle=0$, which means that a substantional number of
eigenstates still have $\tilde{K}^{a}_{\ell}=K^{a}_{\ell}$ with eigenvalues
$\pm 1$.

From the above study on the small system, we found that the modified LIOMs,
$\\{{\tilde{K}}^{a(b)}_{\ell}\\}$, indeed label all energy eigenstates. We
expect that this holds for the system with larger sizes. Also, the interaction
$H_{\rm int}$ mixes only a small number of local cluster-basis eigenstates of
$H_{0}$. This fact implies that the interacting model in the present work
exhibits localization phenomena analogous to Anderson localization. In what
follows, to elucidate it we examine the dynamics of the system.

Numerical demonstration of localization dynamics.— We numerically demonstrate
localization phenomena for larger systems. In what follows, we remove the
‘fictitious’ disorder, $\delta h_{\ell}=0$. The interaction $H_{\rm int}$
affects to the system locally, mixing a few cluster-basis eigenstates, and
also has large kernel space, $H_{\rm int}|\psi_{k}\rangle=0$. From these facts
and the existence of the extensive numbers of LIOMs,
$\tilde{K}^{a(b)}_{\ell}$, we expect localization phenomena in the interacting
system [$H$ in Eq. (3)], especially, the quench dynamics that exhibits non- or
slow-thermalization Rigol2007 .

We observe Loschmidt echo (LE), given by $({\rm
LE})=|\langle\psi(t)|\psi(0)\rangle|^{2}$ where $|\psi(t)\rangle$ is the many-
body wave function at time $t$, and the entanglement entropy (EE), defined as
$S=-\mbox{Tr}[\rho_{A}\ln(\rho_{A})]$ with $A$-subsystem reduced density
matrix $\rho_{A}=\mbox{Tr}_{B}[\rho]$, where $\rho$ is a density matrix of the
entire system and the subsystem is set to $L/2(L/2+1)$-site system for even
(odd) $L$. For the practical calculation, $A$ and $B$ subsystems are set as
shown in Fig. 1. In what follows, time is measured in units $[g/\hbar]$ and we
set a random cluster-basis state
(e.g.,$|1\bar{1}1\bar{1}\bar{1}\cdots\rangle$) as an initial state. In the
quench dynamics, we average over 80 samples for the initial state and disorder
realizations, $\\{\lambda_{\ell}\\}$.

Numerical results are displayed in Fig. 3. The LE remains large finite values
for a long period and the EE remains low-values around $2\ln 2$ clusterEE  as
shown in Figs. 3 (a) and 3 (b) supp1 . These behaviors retain for large $g$’s.
The results of the LE in Fig. 3 (a) indicate that initial-state information is
preserved, which indicates ergodicity breaking and also the behavior of EE
implies all eigenstates of $H$ remains to be low-entangled Supp0 . These
numerical results indicate even for the cluster model with interactions
$H_{\rm int}$, genuine localization dynamics with short-range entanglement
retains, which is similar to Anderson localization.

Figure 3: Dynamics of LE [(a)] and EE [(b)]. We set $L=12$ and averaged over
$60$ disorder realization for $W=2$ and initial random cluster-basis state.
The EE is quite stable, and its value is close to $2\ln 2$ (corresponding to
the value obtained by cut two clusters). We set $\delta h=0$.

Effects of additional interactions: Many-body localization dynamics.— We add
an Ising interaction to the system,
$V_{zz}=\sum^{L-1}_{j=0}v_{zz}\sigma^{z}_{j}\sigma^{z}_{j+1}$, and observe how
the localization nature of the model [Eq. (3)] changes. The interaction,
$V_{zz}$, preserves the $\mathbb{Z}_{2}\times\mathbb{Z}^{T}_{2}$ symmetry, and
therefore, we expect that the SPT order persists even in the presence of
$V_{zz}$. Obviously for finite $v_{zz}$, the operator
$\tilde{K}^{a(b)}_{\ell}$ is no longer exact LIOMs, but some emergent LIOMs
might appear, as suggested in Bahri .

We observe that the level-spacing analysis for the system with finite $v_{zz}$
exhibits localization properties. [See Supp0 ]. Then, we numerically
investigate the quench dynamics of the model with finite $g$ and $v_{zz}$,
where we set the cluster-basis Neel state $|1\bar{1}1\bar{1}1\cdots\rangle$ as
an initial state and observe the LE, EE and a loop order (LO) defined by
$({\rm LO})=\langle\Psi(t)|{\hat{L}}|\Psi(t)\rangle$ with
${\hat{L}}=\prod^{L/2-1}_{\ell=0}K^{a}_{\ell}$ supp2 . The LO diagnoses the
presence of SPT order.

Numerical results are shown in Fig. 4, where we set $g=1$. For small $v_{zz}$,
the value of the LE remains finite [Fig. 4 (a)], the increase of the EE is
much suppressed [Fig. 4 (b)] and also the value of the LO remains finite [Fig.
4 (c)] for a long period. Obviously, these are MBL behavior with SPT order.
For larger $v_{zz}$, localization tendency is weakened, i.e., the increase of
the EE is enhanced with logarithmic growth and saturates with larger values,
approaching to the Page value Page , $(L\log 2-1)/2$. Correspondingly, the
values of the LE and LO also are decreasing, thermalization tendency is
enhanced and the SPT order is fading away. We furthermore investigate the
behavior of a modified LO order described by $\tilde{K}^{a}_{\ell}$, and find
that SPT order again survives for short periods, as the details are shown in
Supp0 .

Conclusion.— We showed that a specific type of interacting cluster model with
disorder exhibits genuine MBL. LIOMs modified from the stabilizer operators
were found for arbitrary strength of interactions and disorder. We showed that
the LIOMs label all energy eigenstate. The locality of the LIOMs and extensive
number of them imply localization phenomena analogous to Anderson
localization. We numerically demonstrated that the quench dynamics of the
system shows non-thermalized dynamics (ergodicity breaking dynamics),
characteristic of the localization phenomena with SPT order. Further, non-
thermalized dynamics is stable against weak interactions. The localization
tendency characterized by the LIOMs [$\tilde{K}^{a(b)}_{\ell}$] survives in
the presence of the Ising interactions, $V_{zz}$.

Finally, we comment that another type of the LIOMs can be defined, which have
not compact but a long-tail support LIOM_long . The study of such modified
LIOM system is an interesting future topic.

Figure 4:  Dynamics of LE [(a)], EE [(b)] and $LO$ [(c)]. We set $L=12$,
$W=2$, $g=1$ and $\delta h=0$ and averaged over $40$ disorder realizations.

The work is supported by JSPS KAKEN-HI Grant Number JP21K13849 (Y.K.).

## References

  * (1) R. Nandkishore, and D. A. Huse, Annual Review of Condensed Matter Physics 6, 15 (2015).
  * (2) D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Rev. Mod. Phys. 91, 021001 (2019).
  * (3) D. M. Basko, I. L. Aleiner, and B. L. Altshuler, Ann. Phys. 321, 1126 (2006).
  * (4) J. H. Bardarson, F. Pollmann, and J. E. Moore, Phys. Rev. Lett. 109, 017202 (2012).
  * (5) M. Rigol, V. Dunjko, V. Yurovsky, and M. Olshanii, Phys. Rev. Lett. 98, 050405 (2007).
  * (6) 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).
  * (7) C. J. Turner, A. A. Michailidis, D. A. Abanin, M. Serbyn, and Z. Papić, Nat. Phys. 14, 745 (2018).
  * (8) S. Choi, C. J. Turner, H. Pichler, W. W. Ho, A. A. Michailidis, Z. Papić, M. Serbyn, M. D. Lukin, and D. A. Abanin, Phys. Rev. Lett. 122, 220603 (2019).
  * (9) W. W. Ho, S. Choi, H. Pichler, and M. D. Lukin, Phys. Rev. Lett. 122, 040603 (2019).
  * (10) T. Iadecola and M. Schecter, Phys. Rev. B 101, 024306 (2020).
  * (11) D. Bluvstein, A. Omran, H. Levine, A. Keesling, G. Semeghini, S. Ebadi, T. T. Wang, A. A. Michailidis, N. Maskara, W. W. Ho, S. Choi, M. Serbyn, M. Greiner, V. Vuletic, and M. D. Lukin, arXiv:2012.12276 [quant-ph] (2020).
  * (12) P. A. McClarty, M. Haque, A. Sen, and J. Richter, Phys. Rev. B 102, 224303 (2020).
  * (13) N. Shibata, N. Yoshioka, and H. Katsura, Phys. Rev. Lett. 124, 180604 (2020).
  * (14) K. Lee, R. Melendrez, A. Pal, and H. J. Changlani, Phys. Rev. B 101, 241111 (2020).
  * (15) Y. Kuno, T. Mizoguchi, Y. Hatsugai, Phys. Rev. B 102, 241115 (R) (2020).
  * (16) M. Serbyn, D. A. Abanin and Z. Papić, Nat. Phys. 17, 675 (2021).
  * (17) J. Jeyaretnam, J. Richter, and A. Pal, Phys. Rev. B 104, 014424 (2021).
  * (18) Y. Kuno, T. Mizoguchi, Y. Hatsugai, Phys. Rev. B 104, 085130 (2021).
  * (19) M. Serbyn, Z. Papić , and D. A. Abanin, Phys. Rev. Lett. 111, 127201 (2013).
  * (20) D.A. Huse, R. Nandkishore, and V. Oganesyan, Phys. Rev. B 90, 174202 (2014).
  * (21) J. Z. Imbrie, J. Sat. Phys 163, 998 (2016).
  * (22) J.Z. Imbrie, V. Ros, and A. Scardicchio, Ann. Phys. (Berlin) 529, 7, 1600278 (2017).
  * (23) H. J. Briegel and R. Raussendorf, Phys. Rev. Lett. 86, 910 (2001).
  * (24) J. K. Pachos and M.B. Plenio, Phys. Rev. Lett. 93, 056402 (2004).
  * (25) W. Son, L. Amico, R. Fazio, A. Hamma, S. Pascazio, and V. Vedral, Europhys. Lett. 95, 50001 (2011).
  * (26) P. Smacchia, L. Amico, P. Facchi, R. Fazio, G. Florio, S. Pascazio, and V. Vedral, Phys. Rev. A 84, 022304 (2011).
  * (27) B. Bauer and C. Nayak, J. Stat. Mech. Theory Exp. P09005 (2013).
  * (28) Y. Bahri, R. Vosk, E. Altman, and A. Vishwanath, Nat. Commun. 6, 7341 (2015).
  * (29) R. Vasseur, A.J. Friedman, S.A. Parameswaran, and A.C. Potter, Phys. Rev. B 93, 134207 (2016).
  * (30) S.A. Parameswaran and R. Vasseur, Reports Prog. Phys. 81, 082501 (2018).
  * (31) K.S.C. Decker, D.M. Kennes, J. Eisert, and C. Karrasch, Phys. Rev. B 101, 014208 (2020).
  * (32) Y. Kuno, Phys. Rev. Research 1, 032026(R) (2019).
  * (33) T. B. Wahl and B. Béri, Phys. Rev. Research 2, 033099 (2020).
  * (34) A. Chan and T. B. Wahl, J. Phys.: Cond. Mat. 32, 305601 (2020).
  * (35) J. Li, A. Chan, and T. B. Wahl, Phys. Rev. B 102, 014205 (2020).
  * (36) J. Kemp, N. Y. Yao, and C. R. Laumann, Phys. Rev. Lett. 125, 200506 (2020).
  * (37) R. Sahay, F. Machado, B. Ye, C. R. Laumann, and N. Y. Yao, Phys. Rev. Lett. 126, 100604 (2021).
  * (38) C. M. Duque, H. Y. Hu, Y. Z. You, V. Khemani, R. Verresen, and R. Vasseur, Phys. Rev. B 103, L100207 (2021).
  * (39) R. Verresen, R. Moessner, and F. Pollmann, Phys. Rev. B 96 101103 (2017).
  * (40) R. Verresen, N. G. Jones, and F. Pollmann, Phys. Rev. Lett. 120 057001 (2018).
  * (41) A. Smith, B. Jobst, A.G. Green, and F. Pollmann, arXiv:1910.05351 (2020).
  * (42) Y. Kuno, T. Orito, and I. Ichinose, New J. Phys. 22, 013032 (2020).
  * (43) T. Orito, Y. Kuno, and I. Ichinose, Phys. Rev. B 101, 224308 (2020).
  * (44) T. Orito, Y. Kuno, and I. Ichinose, Phys. Rev. B 104, 094202 (2021).
  * (45) See supplemental material: (I) Structure of Hilbert space. (II) The dependence of the fictitious disorder (III) System size dependence of the LE. (IV) Level spacing analysis for finite $v_{zz}$. (V) Calculations of a modified LO.
  * (46) We employed the Quspin solver for all numerical calculations: P. Weinberg and M. Bukov, SciPost Phys. 7, 20 (2019); 2, 003 (2017).
  * (47) This value of EE corresponds to the value when the two cluster states are cut.
  * (48) Additional numerical calculations for system size dependence and the dependence of the fictitious disorder are given in Supplemental material Supp0 .
  * (49) The operator LO is related to the string order parameter, which characterizes the bulk SPT order Son ; Smacchia .
  * (50) D. N. Page, Phys. Rev. Lett. 71, 1291 (1993).
  * (51) If in the model of Eq.(3), we change $H_{\rm int}$ to $H^{\prime}_{\rm int}=\sum_{r\neq 0,\ell}V^{r}_{\ell}$, where $V^{r}_{\ell}\equiv ge^{-|r|}N_{\ell+r}(K^{a+}_{\ell}K^{b-}_{\ell}+K^{a-}_{\ell}K^{b+}_{\ell})$, another type of LIOMs can be constructed such as $L^{a}_{\ell}\equiv K^{a}_{\ell}+\frac{1}{2\lambda_{\ell}}\sum_{r\neq 0}V^{r}_{\ell}$ and $L^{b}_{\ell}\equiv K^{b}_{\ell}-\frac{1}{2\lambda_{\ell}}\sum_{r\neq 0}V^{r}_{\ell}$. These LIOMs are not compact but a long-tail support. Hence, this system may induce MBL phenomena, which poses a future work.
  * (52) V. Oganesyan and D.A. Huse, Phys. Rev. B 75, 155111 (2007).

## Supplemental material

###  .1 I. Structure of Hilbert space under ‘fictitious’ disorder

For $L$-site system, where the total dimension of Hilbert space is given by
$N_{D}=2^{L}$, all eigenstates $|\psi_{k}\rangle$ for $H$ of Eq.(3) in the
main text with $g\neq 0$ under weak ‘fictitious’ disorder are classified into
two classes since the interaction acts to only the following four states in
two unit-cells: $|11\rangle_{\ell-1}|1\bar{1}\rangle_{\ell}$,
$|11\rangle_{\ell-1}|\bar{1}1\rangle_{\ell}$,
$|\bar{1}\bar{1}\rangle_{\ell-1}|1\bar{1}\rangle_{\ell}$, and
$|\bar{1}\bar{1}\rangle_{\ell-1}|\bar{1}1\rangle_{\ell}$.

The first class is composed of eigenstates satisfying $H_{\rm
int}|\psi_{k}\rangle=g\sum_{\ell}V_{\ell}|\psi_{k}\rangle=0$, that is, the
eigenstate $|\psi_{k}\rangle$ is a null state for $H_{\rm int}$. On the other
hand, the second class is composed of the ones satisfying $H_{\rm
int}|\psi_{k}\rangle=v_{k}|\psi_{k}\rangle$, where $v_{k}$ is finite real
value. The above classification can be understood by observing how the
interaction acts on the eigenstates of $H_{0}$.

In $N_{D}$ eigenstates of $H$, $2^{L/2+1}$ eigenstates are totally unaffected
by the interaction, $H_{\rm int}$. These eigenstates consist of the following
two categories: (I) in the eigenstate, the state of each unit-cell is given by
$|11\rangle_{\ell}$ or $|\bar{1}\bar{1}\rangle_{\ell}$ (total $2^{L/2}$
eigenstates). (II) in the eigenstate, the state of each unit-cell are given by
$|1\bar{1}\rangle_{\ell}$ or $|\bar{1}1\rangle_{\ell}$ (total $2^{L/2}$
eigenstates). Needless to say, these total $2^{L/2+1}\equiv N^{1}_{D}$
eigenstates are eigenstates of $H_{0}$.

On the other hand, the number of eigenstates of $H$ with finite interaction
energies is obtained by counting the number of eigenstates of $H_{0}$ affected
by the interaction $H_{\rm int}$. This number can be counted as follows: (a)
We consider the eigenstate in which the state of each unit-cell is given by
$|11\rangle_{\ell}$ or $|\bar{1}\bar{1}\rangle_{\ell}$ (total $2^{L/2}$
eigenstates). (b) For each eigenstate, select $k$ unit-cells and change their
states of the unit-cells to $|1\bar{1}\rangle_{\ell}$ or
$|\bar{1}1\rangle_{\ell}$, where $k$ takes from $1$ to $L/2-1$. The total
number of eigenstates obtained in this way is
$\sum^{L/2-1}_{k=1}2^{L/2-k}\binom{L/2}{k}2^{k}=4^{L/2}-2^{L/2+1}\equiv
N^{2}_{D}$, which are affected by the interaction $H_{\rm int}$. These
$N^{2}_{D}$ eigenstates are the second ones. In addition, the total sum of the
eigenstates of the first and second class is $N^{1}_{D}+N^{2}_{D}=N_{D}$.

Furthermore, the interaction $H_{\rm int}$ mixes the eigenstates of the second
class of $H_{0}$, that is, $N^{D}_{2}$ eigenstates of $H_{0}$. However, the
mixing is small and local. That is, the Hamiltonian matrix of $H$ based on the
$N^{2}_{D}$ eigenstates becomes a block matrix with many small blocks. This
implies that the obtained eigenstates of the Hamiltonian matrix $H$ are low-
entangled, where the deviation of the EE of one of the eigenstates of $H_{0}$
is small.

###  .2 II. Effect of fictitious disorder $\delta h_{\ell}$

To understand the structure of the eigenstate of the Hamiltonian $H$ of Eq.
(3) in the main text, we added very small ‘fictitious’ disorder $\delta
h_{\ell}$. We expect that such a small $\delta h_{\ell}$ gives little effect
to the localization nature of the system. As a concrete examination on this
point, we observe the dependence on $\delta h$ of the LE. The numerical
estimation is shown in Fig. 5. The resultant dynamics of the LE is almost
independent of the strength of $\delta h$ where $\delta
h\leq\mathcal{O}(10^{-4}g)$. From this fact, we expect that other physical
observables in the system dynamics are not affected by the ‘fictitious’
disorder $\delta h_{\ell}$.

Figure 5: Dependence on the strength of ‘fictitious’ disorder $\delta
h_{\ell}$ of the LE with $\delta h_{\ell}\in[-\delta h,\delta h]$ (uniform
distributed disorder). We set $L=12$, $W=2$, $g=1$ and $v_{zz}=0$. The initial
state is a random cluster-basis state. We averaged over 80 disorder and random
initial state samples.

###  .3 III. System size dependence

We examine the system size dependence of the dynamics. In particular, we focus
on the LE and calculated the LE for various system sizes. The result is shown
in Fig. 6. From the left panel of Fig. 6, the system size dependence is small.
All LEs remain at some finite value ($\sim 0.5$) for a long period. We expect
that the finite value of the LE also survives for larger system sizes as shown
in the right panel of Fig. 6.

Figure 6: Left panel: System size dependence of the LE for $L=8,10,12$ and
$14$. We set $W=2$, $\delta h=0$, $g=1$, $v_{zz}=0$. We averaged over disorder
and random initial state samples, the numbers of the sample are 100, 80, 60
and 40 for $L=8,10,12$ and $14$, respectively. Right panel: The system size
dependence of the time average of the LE in from $t=0.03$ to $t=10^{3}$.

###  .4 IV. Level spacing analysis for finite $v_{zz}$

To examine the presence of the localization tendency and the integrability of
the system, we employed the level spacing analysis (LSA) for the system with
finite $v_{zz}$ Oganesyan2007 . We diagonalize the Hamiltonian $H+V_{zz}$,
obtain all energy eigenvalues and calculate the level spacing ratio $r_{s}$
defined by $r_{s}=[{\rm min}(\delta^{(s)},\delta^{(s+1)})]/[{\rm
max}(\delta^{(s)},\delta^{(s+1)})]$ for all $s$, where
$\delta^{(s)}=E_{s+1}-E_{s}$ and $\\{E_{s}\\}$ is the set of energy eigenvalue
in ascending order. Then, we calculate the mean level spacing ratio $\langle
r\rangle$, which is obtained by averaging over $r_{s}$ with employing all
energy eigenvalues and also further averaging over disorder realizations for
$\\{\lambda_{\ell}\\}$. The result for various system sizes and $v_{zz}$ is
shown in Fig. 7. For small $v_{zz}$, the remnant of degeneracy causes the mean
value of the level spacing ratio to be smaller than that of the Poisson
distribution, $\sim 0.384$. However, for larger $v_{zz}$, the value is getting
slightly larger than the Poisson distribution but stays near the value of the
Poisson distribution. This indicates that the system is in a localized phase.

Figure 7:  Mean level spacing ratio $\langle r\rangle$ for the system with
$L=8-14$, $W=2$, $g=1$ and $\delta h_{\ell}=0$. For small $v_{zz}$, $\langle
r\rangle$ is close to the value of Possion distribution. The numbers of the
disorder sample are $100$, $60$, $40$ and $30$ for $L=8$, $10$, $12$ and $14$,
respectively. The red and blue dashed line are $\langle r\rangle=0.5295$ and
$0.386$, which are the values for the Wigner-Dyson and Poisson distributions,
respectively. Figure 8: The dynamical behavior of the MLO for various
$v_{zz}$. We set $L=12$, $W=2$, $g=1$ and $\delta h_{\ell}=0$, and averaged
over 40 disorder realizations for $\\{\lambda_{\ell}\\}$.

###  .5 V. Calculation of an modified loop order

In the study of quench dynamics with finite $g$ and $v_{zz}$, we further
investigate the behavior of a modified loop order defined by $({\rm
MLO})=\langle\Psi(t)|\biggl{[}\prod^{L/2-1}_{\ell=0}\tilde{K}^{a}_{\ell}\biggr{]}|\Psi(t)\rangle$.
Note that the operator depends on the set of disorder. We set the cluster-
basis Neel state $|1\bar{1}1\bar{1}1\cdots\rangle$ as an initial state, which
has $({\rm MLO})=1(-1)$ for $L/2$ is even (odd). The dynamics of the MLO for
various $v_{zz}$ is shown in Fig. 8. For short-time, the values of the MLO
remains the initial values. The initial value tends to decay after a long
period of the time evolution, with large swings. This is due to the fact that
state mixing with large LIOM eigenvalues occurs in the process of long time
evolution.


# Prototype Cherenkov Detector Characterization for Muon Tomography
Applications

T. Avgitas<EMAIL_ADDRESS>J.-C. Ianigro<EMAIL_ADDRESS>J.
Marteau<EMAIL_ADDRESS>
###### Abstract

Muography is an innovative imaging technique using naturally produced
elementary particles – atmospheric muons – like the X-rays of medical imaging.
The modification of the particles flux – by scattering or absorption –,
reflects the contrasts in density within the medium and therefore offers the
possibility for an image of the crossed volumes. The imaging process is based
on the tracking of the particles which accounts for the absorption or the
scattering of the muons trajectories. Neither the energy nor the identity of
the particles (the so-called PID) is exploited since this information
traditionally relies on the use of calorimeters and/or high intensity magnetic
fields. Both these techniques hinder detector portability which in the case of
muography is important and this renders them impractical for its purpose. In
this paper we characterize the performance of a simple and small water
Cherenkov detector capable on the one hand of providing some insights on
energy and PID and on the other hand of improving the background rejection for
a muon telescope. We tested a prototype of such water Cherenkov detector in
combination with two small muon hodoscopes. Both systems are using the same
opto-electronics chain – optical fibers and pixellized photosensors – and the
same data acquisition (DAQ) readout system which ensures an easy integration
and implementation within presently running systems. This article presents the
test setup, the detector response to cosmic muons and its performance
evaluation against a basic simulation of its geometry and detection principle.

###### keywords:

muography , Cerenkov , tracking , PID , background rejection

††journal: Nuclear Physics B

[affil]organization=IP2I, Univ Lyon, Univ Claude Bernard Lyon 1, CNRS/IN2P3,
IP2I Lyon, UMR 5822, addressline=, postcode=F-69622, city=Villeurbanne,
country=France

##  1 Introduction

Muon imaging or muography has emerged as a powerful non-invasive method to
complement standard tools in Earth Sciences and is nowadays applied to a
growing number of fields such as industrial controls, homeland security, civil
engineering. This technique relies on the detection of modifications -
absorption or scattering - in the atmospheric muon flux when these particles
cross a target. Atmospheric muons are secondary products of primary cosmic-
rays, essentially protons and helium nuclei expelled by stars, interacting
with nuclei encountered on the top of the atmosphere. The rather low
interaction cross-section of muons with matter ensures that most of them reach
the Earth’s ground level and that furthermore they may significantly penetrate
large and dense structures. As suggested originally by Alvarez in 1970 for the
Chephren pyramid, this property may be exploited to perform density contrasts
analysis of the interior of the target like X-rays do in medical imaging.
Details on the basic muography features may be found in (Marteau et al., 2012)
and references therein. Essentially muography relies on the tracking
capabilities of portable detectors that are positioned behind (or on either
side) of the structure under study. Many detection technologies are used so
far : scintillators (Lesparre et al., 2012), emulsions (Tioukov et al, 2019) ,
resistive plate chambers (Ambrosino et al., 2015), micro-megas (Bouteille et
al., 2016) etc. They have various performances in spatial and angular
resolutions which are in turn used to assess the absorption or the scattering
efficiencies and then, after the resolution of an inverse problem, the density
distributions in the target. But the tracking efficiency is not sufficient in
principle when one aims to solve an inverse problem in muography, which is by
essence constraint by the physics of muon-matter interaction. For instance,
inverting data in scattering muography involves the momentum $p$ of the
particle which enters in the multiple Coulomb scattering formula :

$\theta_{0}=\frac{13.6}{\beta
cp}z\sqrt{\frac{x}{X_{0}}}\left[1+0.0038\ln\left(\frac{x}{X_{0}}\right)\right].$

Inferring the particle momentum or energy is therefore a way to constrain the
inverse problem algorithm to select the best solution among all the possible
ones. Since the measurement of energy/momentum of high energy particles
requires detection capabilities well beyond the real possibilities of field
equipment powered through solar panels for example. Exploiting the Cerenkov
effect has been suggested already in the past (Peña et al., 2021) on the model
of the AUGER experiment Abreu et al (2021). The idea is to use water-based
Cerenkov that may also provide interesting background rejection capabilities.
In the MUTE experiment for example, a large water cuve is placed downstream
from 2 scintillator planes. The typical dimensions of the cuve are of the
ordre of 1m. The Cerenkov induced photons are collected by a large
photomutiplier tube (PMT) while the scintillator planes are readout with
SiPMs. In this article we present first results obtained with a prototype of
water Cerenkov detector where the light collection is performed by a bundle of
WLS optical fibers connected via a mechanical cookie to a multi-anode PMT.
This design offers multiple benefits among which the use of a compact readout
system with no large size Photomultiplier Tube (PMT) such as the ones used in
the AUGER experiment for instance. It allows also to re-use the same readout
system as the one for the scintillators planes performing the tracking, with
plastic scintillators planes of various transverse segmentation and active
surface. This readout uses a custom, low power consumption, auto-triggerable
Ethernet capable system described for instance in (Marteau et al., 2012).
Details are also given below. The paper describes the setup and the results
obtained. It is organized as follows. The second and third sections describe
the Cerenkov prototype and test setup. The fourth section details the track
selection procedure while the fifth one gives the general outlines of the
simulation tool developed for the data/MC comparison. Last section gives the
main results and comment on the reasonable agreements achieved and opens on
the perspectives for such detector design.

##  2 Cherenkov detector description

For test purposes we designed the simplest possible tank for the Cherenkov
detector prototype that consists of a cylinder of grey PVC with 1 cm thickness
and 15.7 cm radius. The internal surface is covered with Tyvek to increase
reflectivity, the water level inside stands at 1 m. A Hamamatsu H8804-300
Multianode Photomultiplier Tube (MaPMT) is used for the photon detection. The
MaPMT window is optically coupled to 64 optic fibers (Kuraray-Y11) the length
of which runs parallel to the cylinder’s axis (fig. 1 \- left), their purpose
is to increase the light collection efficiency by delivering the cherenkov
photons to the MaPMT window.

Figure 1: (left) The Cherenkov detector light collection system of fibers. The
experimental setup (center) and its simulation counter part (right) for
testing the DAQ and the Cherenkov response to atmospheric muons.

The submerged part of the fibers is 80 cm in length, it has a rope-like
structure with the fibers tangling with each other, something that should give
rise to shadowing effects for the inner-most fibers of the bundle hindering
collection efficiency and should be optimized in the future. This attribute is
not a problem here since the purpose of this work was to test the feasibility
of the DAQ, test how it operates in time coincidence with scintillation
detectors and get insights on the capabilities, the potential drawbacks and
guide future developments of the Cherenkov detector itself.

### Characterization setup

Two muon hodoscopes are placed on either side of the cherenkov detector, one
top and one below (fig. 1 \- center/right). The top hodoscope is enclosed in a
metal box to facilitate the alignment of the two detection planes. Each
hodoscope consists of two detection planes that monitor the traversing
particles’ positions ($x$, $y$) with scintillator strips. The $z$ coordinate
of the particle hits is provided by the planes’ height positions since all
four of them run parallel to the ground.

The top two planes consist of 16 scintillator strips per $x$ and $y$ direction
with strip dimensions 2.5 cm $\times$ 40 cm and a distance of $\Delta z$ = 34
cm between planes while the bottom two planes are built with 8 scintillator
strips per direction with 5 cm $\times$ 40 cm dimensions and set at a distance
of $\Delta z$ = 49.5 cm between them. The DAQ setup and the readout
electronics are similar to those used for the Diaphane telescopes (Lesparre et
al., 2012), specifically for every detector, scintillator or Cherenkov, the
photons are guided by WLS fibers to the window of a MaPMT. One MaPMT monitors
the 64 scintillator strips for the top hodoscope and another one the 32 strips
for the bottom hodoscope.

The goal of this work is to study the Cherenkov detector response to muon
tracks crossing the scintillator surfaces, to develop a modelization of the
detector and test in the future its potential energy and particle
identification capabilities. To investigate this we performed two runs, one
with lead bricks of 40 cm total height, between the bottom hodoscope planes
and one without lead that we used as a reference. The duration of the first
measurement (reference) was 453 hours and the second run, with the lead
bricks, lasted 553 hours.

### Triggering and track selection

The top hodoscope operates with an “and” trigger-0 condition between planes
which means it triggers only when both planes register signal. The bottom
hodoscope’s trigger-0 condition is “or” which means that triggers when at
least one out of the two planes registers signal. The trigger-0 condition for
the Cherenkov detector is a simple voltage threshold applied by the
electronics to the MaPMT signal. The trigger-1 condition that triggers the
data registration for the entire experimental setup is set in “or” mode, which
means that at least two trigger-0 conditions need to be satisfied. All the
data discussed onward are only trigger-1 events. The reason behind this
triggering schema is that we are interested in differentiating between
particles that reach the rear bottom plane and those which stop in between
planes for the bottom hodoscope (particle energy threshold detection).

The muons we are interested in are those that cross the entire volume of water
inside the Cherenkov detector from the top water level to its bottom. We call
these tracks “contained” to denote that they remain within the cylinder and
they don’t escape from the sides. The containment is determined from the
interpolation of the reconstructed particle tracks to the height of the water
level and the reconstructed track position on the detection plane 1 (fig 1 \-
right) placed just below the cylinder.

### Event selection

A particle that crosses a detection plane should appear as a “hit” with ($x$,
$y$) coordinates. This is not the case because at the level of the MaPMT there
is crosstalk between channels and a particle crossing can trigger multiple
channels. To mitigate this effect for the top MaPMT we apply a charge cut for
all channels at the “valley” where the tail of the noise charge meets the rise
of the photo-electron distribution. For the bottom MaPMT the DAQ threshold is
high enough to cut these low charges. We treat the four detection planes of as
unique hodoscope. The event selection process gives rise to three event
categories:

  * 1.

quads, where all detection planes had at least one hit.

  * 2.

triples1, where the three top planes had hits (Det. 1, 2, 3 - fig. 1)

  * 3.

triples2, where the top two planes and the rear bottom one had hits (Det. 0,
2, 3 - fig. 1)

### Track selection

The track reconstruction for each event is done by performing a separate
linear regression fit of the hits’ projections on the $xz$ and the $yz$
planes. Events have different multiplicity of hits per hodoscope plane which
lead to multiple track reconstructions for the same event. We define as a
track the combination of hits that has the minimum chi-squared value for each
projection plane. The chi-squared histograms for both directions ($x$ vs $z$
and $y$ vs $z$) are not all the products of tracks and they include also large
values. We selected chi-squared cut value at the point where the distribution
for the track populations change slope (in semi-log scale) which is around
7.81 for quads and at 5.99 for both triples populations.

The reconstructed tracks sometimes happen to pass outside of the muon
hodoscope borders. This in a small portion of cases is caused by inaccuracies
in the exact positioning of the detectors which affect the reconstruction of
the track. The most important effect though comes from the tracks of triples1
events. Extrapolating these tracks to the level of the last detector shows
that the vast majority of these are actually tracks that are not detected
because of the dimensions of the last detection plane. We mention this because
there are also tracks that fall inside the last detection plane when
extrapolated but are not detected due to the detection efficiency. We apply
cuts for the reconstructed tracks position at the height of each detection
plane from 0 to 40 cm in order to minimize this geometric effect.

The detection inefficiency is $\sim$10% for each detection plane of the bottom
hodoscope. This means that for the both measurements, with or without lead,
triples2 events are due to the detection inefficiency of Det1 and for this
reason we incorporated them into the quads sample. The same holds true for the
triples1 measurements done with lead but it is impossible to separate between
the part of the sample of triples1 that is due to the inefficiency of the Det0
detector and that caused by muons without enough energy to cross the lead
bricks. We keep the triples1 tracks sample separate from the quads sample also
for the measurements performed without lead to highlight the difference (see
sec. 4).

##  3 Simulation

### Detector geometry

We developed a simplified simulation based on GEANT4 libraries (Agostinelli et
al. (2003), Allison et al. (2006), Allison et al. (2016)) to evaluate the
expected outcome for our characterization setup. The simulated muon flux
assumes perfect alignment of the setup constituents around the $z$ axis,
isotropy, and a complete shadowing of our detectors by the building floors on
top. The first two assumptions have been already discussed to hold true only
in approximation, the third one is true for the small zenithal angles
corresponding to the contained tracks we use. The simulation is very
simplistic but since the goal of this study was to test the DAQ setup and
evaluate its performance for a Cherenkov detector we don’t consider this a
drawback.

The simulated counterpart of the Cherenkov detector respects the geometric
characteristics but fails to incorporate the optical properties of the
materials. There is no reflection, refraction nor absorption for the photons
traversing the water or coming in contact with the Tyvek surfaces. The same
holds true for the optical fibers. The only thing taken into account is the
wavelengths of the cherenkov photons and every photon that touches the optic
fibers is considered detected provided that its wavelength falls within the
limits of the absorption spectrum of the fibers (fig. 2(a)). The distribution
of the fibers is ideal, nothing like the true tangled geometry inside the
actual detector, it respects the interaxial distance between fibers dictated
by the mounting cookie throughout their entire length of 80 cm. The core of
the optical fibers is made of polystyrene. We modeled the polystyrene
refractive index dispersion
($n^{2}-1=\frac{1.4435\lambda^{2}}{\lambda^{2}-0.02016}$, Sultanova et al.
(2009), fig. 2(b)) so that we retrieve also the Cherenkov photon contribution
from particles crossing the fibers themselves. This effect should be more
pronounced for the actual entangled geometry of the fibers.

(a) Y11 - WLS absorption spectrum

(b) Polystyrene refractive index

Figure 2: (left) Absorption and Emission Spectra for the WLS fiber used for
the light collection of the Cherenkov detector
(https://www.kuraray.com/uploads/5a717515df6f5/PR0150_psf01.pdf). (right)
polystyrene refractive index n (from https://refractiveindex.info/).

The muon flux simulation starts by producing Geantinos uniformly on the top
plane of the bottom hodoscope with an isotropic azimuthal distribution. The
zenithal angle is considered smaller than $None$ (based on the setup geometry)
and the minimum energy to cross the building is 0.597 GeV. For the muon
production we use the analytic formula found in (Shukla et Sankrith, 2018).
The Geantinos propagate upwards toward the roof of the building, they first
need to cross a surface a little larger than that of detector 3 before finally
reaching the area of the roof where they are “killed” and the information
about the position, the direction and their energy is stored. The muon
production incorporates this information to propagate them back towards the
setup at the basement. For the analysis we use tracks that triggered both
planes of the top hodoscope and at least one detection plane of the bottom
hodoscope. For the simulation we use a trigger condition based on this and we
include the Cherenkov detector response only for these events.

### Detector response

The simulation shows that the Cherenkov detector’s signal has a sudden peak
followed by a very long tail. The peak can be described by a Gaussian and the
largest proportion of the tail by an exponential but this needs to be studied
further. Here we use a convolution of a Landau with a Gaussian distribution
(a.k.a. Langau111https://root.cern.ch/root/html404/examples/langaus.C.html)
since this is the function we use also for the experimental data and we want
to be consistent with our comparisons, in the future we will need to come up
with a model for this. The first remark is that the 40 cm of lead have a small
effect on the muon absorption, setting a minimum energy threshold for quads at
$\sim$0.7 GeV, $\sim$500 MeV (from lead) + $\sim$200 MeV from water (Groom et
al., 2001). The expected detector response for these two energy intervals
(below and above 0.7 GeV) is shown in fig. 3 in number of detected photo-
electrons under the assumption of perfect efficiency.

Figure 3: Cherenkov detector simulated response to two different muon energy
intervals above (left) and below (right) 0.7 GeV.

We can safely assume from this that the differences for the Cherenkov response
between contained quads and triples1 in the lead configuration will also come
from differences in the muons energy and they will not be the product of just
different angular distributions for the two populations. It also shows that
there is a difference of 1000 photons for the mean values of the two
histograms which becomes a more modest difference of 300 photons for the fit
of the peaks. This may provide a weak insight on the energy of the particles
and needs to be studied properly for the final detector geometry.

##  4 Data/MC comparison

### Cherenkov detector response

In table 1 we codify an overview of the results for the different measurements
and the simulation. The key quantity we use is the most probable value
position for the Langau. The percent column is the percentage between the
number of events for a certain type and the total events detected. On the
simulation level it gives some insight to the low energy particles effect on
the Cherenkov detector response. On the experimental side the no-lead
measurements show an inefficiency of the bottom hodoscope detector plane
(Det1) at around 10%. The triples2 for the experiment (already discussed) are
the result of the Det0 efficiency and for both measurements sets are
incorporated into the quads sample. The simulation shows clearly that the
triples1 correspond to the low energy part of the muon spectrum. We see that
the Langau peak is lower for these low energy muons and higher for high energy
muons (quads) in both simulation and experiment.

Detector Response for Contained Tracks - Experimental
---
| No Lead | Lead
Type |  Events |  Percent |  Peak Position |  Events |  Percent |  Peak Position
quads |  44976 |  88$\%$ |  29.95 $\pm$ 0.05 pe |  43405 |  71$\%$ |  29.57 $\pm$ 0.05 pe
triples1 |  6290 |  12$\%$ |  29.06 $\pm$ 0.05 pe |  17929 |  29$\%$ |  26.62 $\pm$ 0.08 pe
total |  51266 |  |  61334 |
Detector Response for Contained Tracks - Simulation
| No Lead | Lead
Type |  Events |  Percent |  Peak Position |  Events |  Percent |  Peak Position
quads |  951439 |  99.88$\%$ |  2939 $\pm$ 5 pe |  929305 |  97$\%$ |  3050 $\pm$ 5 pe
triples1 |  1145 |  0.12$\%$ |  2710 $\pm$ 80 pe |  25746 |  3$\%$ |  2800 $\pm$ 22 pe
total |  952584 |  |  955051 |
Table 1: Cherenkov detector response Langau peak position for different types
of tracks - Simulation & Experiment.

### Angular distribution of tracks

We provide the zenithal angular distribution ($\theta=0\deg$ for vertical
muons) for the muon tracks without lead in two representations (fig. 4). We
see in both plots that the experimental data are skewed towards higher values,
an effect that is more pronounced in the right plot. We attribute this
divergence of the experimental data to the distribution of matter for the
building above which is very difficult to implement properly for the Monte-
Carlo simulation. The overall normalization though is correct.

Figure 4: Contained Tracks - No Lead - Simulation Vs Experimental Data -
Angular Distributions.
Simulation (blue line) and experiment (red crosses) for the zenithal
distribution (left) and the azimuthal distribution (right).

### Low energy events

In the experimental setup with lead between the bottom hodoscope’s detectors
the response of the cherenkov detector to quads (& triples2) corresponds to
muons with higher energies (>0.7 GeV) while the triples1 correspond to the
lower part of the spectrum. The ratio of the cherenkov response to triples1
and quads shows (fig. 5) that the detector differentiates for the low energy
part of the spectrum and the high energy part of the spectrum at 26 - 30 pes,
where the charge ratio becomes $\sim$1\. The probability that the cherenkov
response to low energy muons is lower than 26 pes is 1.5 to 2 times higher
than for the high energy tracks to give signal at that region. This result is
consistent with the simulation expectation.

Figure 5: Contained Tracks - No Lead - Simulation Vs Experimental Data -
Angular Distributions.
Simulation (blue line) and experiment (red crosses) for the zenithal
distribution (left) and the azimuthal distribution (right).

We see a similar behavior but reversed and less exaggerated for the charge
response region above 55 pes. Here the probability (on a bin by bin basis)
that the lower energy muons would give signal at this region is 20% less than
for the high energy ones. The behavior of the experimental results for the
cherenkov response region is also consistent with the expectations from the
Monte-Carlo.

Finally there is a region, between 30 pes and 50 pes, for which the agreement
between simulation and experiment is only qualitatively similar. The
simulation shows that there is a very small probability of differentiating
between low energy and high energy events in contrast to the detector that
seems to give up to 2 times more chance of the high energy particles to give
charge in this region. This is something that we need to investigate further
for the future iteration of the cherenkov detector design.

##  5 Conclusions

This article presents the design of a portable, compact, robust and low power
water-Cerenkov detector to be used in field muography experiments. The results
obtained show that Cerenkov photons may be collected by a bundle of WLS fibers
and amplified by the same MaPMT than those used for the scintillator readout.
This simple particular design offers many possibilities for the fibers layout
in the water tank and for the integration of the Cerenkov sub-detector in any
kind of scintillator tracker. The use of the Cerenkov information gives access
to inputs on the incident particle energy, as we demonstrated with a simple
high-energy pass filter made of lead bricks in our setup. The present study
allowed to validate a design and its Monte-Carlo model. These inputs may be of
primary importance to better constrain the reconstruction algorithms of the
muography imaging especially in scattering mode. Next steps are the production
of an optimized detector where the layout of the fibers is upgraded to have a
maximal Cerenkov light collection efficiency for downward propagating
particles (ie the muon crossing the target) and a detailed study of the
muon/electron separation power. The newly designed Cerenkov detector will be
implemented on running muon stations where it will replace the passive lead
shielding used to remove low-energy background at present. This implementation
fits with the present mechanical frame and with the already running DAQ
system.

## Acknowledgements

This work has been carried out with the support of the LabEx LIO of the
Université Claude Bernard Lyon 1, which has been created as part of a program
for future development (reference number ANR-10-LABX-066) initiated by the
French government and overseen by the Agence Nationale de la Recherche (ANR).
This study is part of the ANR-19-CE05-0033 MEGAMU project.

## References

  * Marteau et al. (2012) Marteau, J. and Gibert, D. and Lesparre, N. and Nicollin, F. and Noli, P. and Giacoppo, F., Muons tomography applied to geosciences and volcanology, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment, Science Direct, 2012.
  * Lesparre et al. (2012) Lesparre, N., Marteau, J., Déclais, Y., Gibert, D., Carlus, B., Nicollin, F., and Kergosien, B., Design and operation of a field telescope for cosmic ray geophysical tomography, Geosci. Instrum. Method. Data Syst., 1, 33–42, https://doi.org/10.5194/gi-1-33-2012, 2012.
  * Tioukov et al (2019) Tioukov, Valeri and Alexandrov, Andrey and Bozza, Cristiano and Consiglio, Lucia and D’Ambrosio, Nicola and De Lellis, Giovanni and De Sio, Chiara and Giudicepietro, Flora and Macedonio, Giovanni and Miyamoto, Seigo and others, First muography of Stromboli volcano, Scientific reports, 9, 1, 6695, 2019, Nature Publishing Group UK London.
  * Ambrosino et al. (2015) Ambrosino, F and Anastasio, A and Bross, A and Béné, S and Boivin, P and Bonechi, L and Cârloganu, C and Ciaranfi, R and Cimmino, L and Combaret, Ch and others, Joint measurement of the atmospheric muon flux through the Puy de Dôme volcano with plastic scintillators and Resistive Plate Chambers detectors, Journal of Geophysical Research: Solid Earth, 120, 11, 7290–7307, 2015, Wiley Online Library.
  * Bouteille et al. (2016) Bouteille, S. et al., A Micromegas-based telescope for muon tomography: The WatTo experiment. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 834, 223-228 (2016).
  * Peña et al. (2021) Jesús Peña-Rodríguez, Luis A. Núñez, Hernán Asorey, Characterization of the muography background using the Muon Telescope (MuTe).
  * Abreu et al (2021) P. Abreu _et al._ (Pierre Auger Collaboration), http://dx.doi.org/10.1088/1748-0221/7/10/P10011, J. Instrum. 7, P10011 (2012)
  * Agostinelli et al. (2003) S. Agostinelli et al., Geant4 - A Simulation Toolkit, Nucl. Instrum. Meth. A 506 (2003) 250-303
  * Allison et al. (2006) J. Allison et al., Geant4 Developments and Applications, IEEE Trans. Nucl. Sci. 53 (2006) 270-278
  * Allison et al. (2016) J. Allison et al., Recent Developments in Geant4, Nucl. Instrum. Meth. A 835 (2016) 186-225
  * Sultanova et al. (2009) N. Sultanova, S. Kasarova and I. Nikolov. Dispersion properties of optical polymers, Acta Physica Polonica A 116, 585-587 (2009) (fit of the experimental data with the Sellmeier dispersion formula: Mikhail Polyanskiy)
  * Shukla et Sankrith (2018) P. Shukla and S. Sankrith, Energy and angular distributions of atmospheric muons at the Earth, International Journal of Modern Physics AVol. 33, No. 30, 1850175 (2018)
  * Groom et al. (2001) D. E. Groom, N.V. Mokhov, S. I. Striganov, Muon stopping power and range tables 10 MeV-100 TeV, Atomic Data and Nuclear Data Tables, Vol. 76, No. 2, July 2001