# Higher-order nonclassicalities in hybrid coherent states

Sandip Kumar Giri1, Biswajit Sen2<EMAIL_ADDRESS>
###### Abstract

In the development of quantum technologies, nonclassical states have been
playing a pivotal role, as quantum advantage cannot be obtained without
appropriate utilization of nonclassicality. In the present work, we consider a
hybrid coherent state (HCS), which is a coherent superposition of the single-
photon-added coherent (SPAC) state and a coherent state (CS). Here, we report
higher-order nonclassical properties of HCS with a specific focus on higher-
order squeezing and higher-order antibunching. It’s shown that HCS is
experimentally realizable, and this engineered quantum state can be used to
produce quantum states with desired higher-order nonclassical properties.

1Department of Physics, Panskura Banamali College, Panskura-721152, India

2Department of Physics, Vidyasagar Teachers’ Training College,
Midnapore-721101, India

##  1 Introduction

The field of quantum state engineering has been considerably advanced in the
recent past, and several engineered quantum states have been produced [1, 2].
Further, applications of the nonclassical engineered states have been reported
in quantum teleportation, quantum communication, and quantum information
theory [3]. Over the last few decades, various schemes have been proposed to
develop different kinds of quantum states, including squeezed state, Fock
state, etc. Furthermore, some new kinds of states that are generated by adding
or subtracting a photon with a coherent state have been found to show very
interesting properties. For example, we may mention that the photon-added
coherent state was first introduced by Agarwal and Tara [4] in 1991. It is a
state intermediate between the Fock state and the coherent state, and it shows
various nonclassical features [4, 5] like photon antibunching, squeezing, etc.
Now, the single-photon-added coherent state, indicated by the symbol
$\left|\alpha,1\right\rangle$ is produced when a photon is added to a coherent
state. Such a state was experimentally produced by Zavatta et al. [6] using
basic optical processes such as beam splitting, frequency down-conversion [7],
and homodyne detection. Due to the distinctive characteristics of the SPAC
state, it has a variety of uses in quantum optics and quantum information
technology, such as quantum key distribution [8, 9], quantum communication
[10], quantum sensing [11], and quantum metrology [12]. The SPAC state serves
as a versatile test-bed for studying the strange nature of quantum mechanics.
It shows the nonclassical features, such as quadrature squeezing [13],
negativity of the Wigner function [14], sub-Poissonian photon statistics [4],
etc. Recently, Yususf Tarek et al. [15] proposed a new kind of state termed
HCS, which is generated by adding a coherent state to a single-photon-added
coherent state. The HCS has several advantages over the SPAC state [15]. The
nonclassical features of HCS, such as photon statistics, quadrature squeezing,
and amplitude squared squeezing, have been extensively studied [15], although
the higher-order nonclassicalities of HCS have not been studied yet. It’s
known that the higher-order nonclassicalities of engineered coherent states
often exhibit more pronounced nonclassical effects than their lower-order
counterparts, and they can be used for various quantum information processing
tasks, such as quantum computation, quantum communication, quantum
cryptography, etc. [16]. For example, higher-order squeezing is shown to
enhance different aspects of quantum computing technology, such as quantum
communication and imaging [17, 18, 19]. Keeping these in mind, we aim to study
the higher-order nonclassical effects (e.g., higher-order Hong-Mandel type
squeezing, higher-order antibunching) in HCS.

The rest of the paper is organized as follows: As this paper aims to study
some nonlclassical properties of HCS, Section 2 is dedicated to the
description of HCS. Subsequently, in Section 3, a method for experimental
generation of HCS is described with an appropriate schematic diagram of the
possible setup to establish that the results of the present work are
experimentally realizable with the existing technology. In Section 4, it’s
shown that HCS exhibits higher-order squeezing, and the amount of maximum
squeezing increases with the order of squeezing. A similar nature is observed
for higher-order antibunching and the same is described in Section 5. Finally,
the paper is concluded in Section 6.

##  2 Hybrid coherent state

An interesting engineered quantum state is HCS [15], which is a superposition
between a closest to classical state, called a coherent state, and a
nonclassical (SPAC) state where the amplitudes of the CS and the CS on which a
photon is added (equivalently, the creation operator $a^{\dagger}$ is applied)
are the same. Such a superposition can be written as

$\begin{array}[]{lcl}\left|\varPsi\right\rangle&=&c_{1}\left|\alpha\right\rangle+c_{2}e^{i\varphi}\left|\alpha,1\right\rangle,\end{array}$
(1)

where $\left|\alpha\right\rangle$ and $\left|\alpha,1\right\rangle$ represent
CS and SPAC state, respectively, and $c_{i\in\left\\{1,2\right\\}}$are their
amplitudes, respectively. The CS and the SPAC states are also separated by a
phase angle $\varphi$. For $c_{1}=0$ and $c_{2}=0$ the state
$\left|\varPsi\right\rangle$ reduces to SPAC state and CS, respectively. Thus,
it can be treated as an intermediate state. The SPAC state is a well-studied
nonclassical state, and it’s easy to visualize its nonclassical nature as the
addition of a photon leads to the creation of a hole in the photon number
distribution, which is a witness of nonclassicality [20, 21]. The question is,
whether the nonclassicality present in the SPAC state remains preserved in its
superposition with CS? We will address this question in this paper. However,
before we do so, we would like to note that the combined properties of CS and
the SPAC state in HCS would depend on the values of the superposition
parameters $c_{1}$ and $c_{2}$. Also, the SPAC state reduces to a single-
photon state for $\left|\alpha\right|\rightarrow 0$ and a coherent state for
$\left|\alpha\right|\rightarrow\infty$. The coherent state
$\left|\alpha\right\rangle$ is the eigenstate of the annihilation operator $a$
with eigen value $\alpha$, i.e.,
$a\left|\alpha\right\rangle=\alpha\left|\alpha\right\rangle.$ The complex
eigen value $\alpha=\left|\alpha\right|e^{i\zeta}$, where $\zeta$ is a phase
angle and a real number. The coherent state in the Fock basis
$\left|n\right\rangle$ can be represented as

$\begin{array}[]{lcl}\left|\alpha\right\rangle&=&e^{-\frac{\left|\alpha\right|^{2}}{2}}\sum_{n=0}^{\infty}\frac{\alpha^{n}}{\sqrt{n!}}\left|n\right\rangle\\\
&=&e^{\left(\alpha a^{\dagger}-\alpha^{*}a\right)}\left|0\right\rangle\\\
&=&\mathcal{D\left(\alpha\right)}\left|0\right\rangle,\end{array}$ (2)

where $\mathcal{D\left(\alpha\right)}$ is the displacement operator.

##  3 Possible experimental generation of a hybrid coherent state

The HCS can be prepared through the coupling of a single-photon Mach-Zehnder
interferometer with a nonlinear optical Kerr medium [15, 22]. A single-photon
state $\left|1\right\rangle$ and a vacuum state $\left|0\right\rangle$ are
injected into the two input ports of the BS1, which is a 50:50 beam splitter
as shown in Fig. 1. The resulting state at the output ports $b$ and $c$ of BS1
will be

$\left|Out\right\rangle_{BS1}=\frac{1}{\sqrt{2}}\left(\left|10\right\rangle_{bc}+i\left|01\right\rangle_{bc}\right),$
(3)

Figure 1: The scheme for the generation of HCS

where we have used the obvious notation
$\left|1\right\rangle_{b}\left|0\right\rangle_{c}\equiv\left|10\right\rangle_{bc}$.
The state in arm $b$ picks up a phase shift $\theta$. So, after phase
shifting, we have the state

$\left|Out\right\rangle_{BS1,\theta}=\frac{1}{\sqrt{2}}\left(e^{i\theta}\left|10\right\rangle_{bc}+i\left|01\right\rangle_{bc}\right).$
(4)

A coherent state $\left|\alpha\right\rangle$ is present in arm $a$. So, before
the cross-Kerr interaction, we have the input state as

$\left|In\right\rangle_{ck}=\left|Out\right\rangle_{BS1,\theta}\left|\alpha\right\rangle_{a}.$
(5)

Now, the state in arm $b$ is coupled with a coherent beam
$\left|\alpha\right\rangle$ in arm $a$ through a nonlinear cross-Kerr medium,
which induces a cross-Kerr phase modulation (XPM). The related interaction
Hamiltonian is given by $H_{ck}=\hbar Ka^{\dagger}ab^{\dagger}b$, where $K$ is
proportional to $\chi^{(3)}$, the third-order nonlinearity of the Kerr medium.
The action of the Kerr medium yields the output state as

$\left|Out\right\rangle_{ck}=e^{-\frac{iH_{ck}t}{\hbar}}\left|In\right\rangle_{ck}=e^{-i\varPhi
a^{\dagger}ab^{\dagger}b}\left|In\right\rangle_{ck},$ (6)

where $\varPhi=Kt$ is the phase shift of XPM in the interaction time $t$. For
weak XPM, the phase shift angle $\varPhi<<1$. Under such condition, we can
approximate $e^{-i\varPhi a^{\dagger}ab^{\dagger}b}\approx\left(1-i\Phi
a^{\dagger}ab^{\dagger}b\right)$. So, under the weak XPM, we have the output
state as

$\left|Out\right\rangle_{ck}=\frac{1}{\sqrt{2}}\left[\left(1-i\Phi\alpha
a^{\dagger}\right)e^{i\theta}\left|10\right\rangle_{bc}\left|\alpha\right\rangle_{a}+i\left|01\right\rangle_{bc}\left|\alpha\right\rangle_{a}\right].$
(7)

Now, the BS2 is a variable beam splitter, which causes the transformations

$\begin{array}[]{rcl}\left|10\right\rangle_{bc}&\stackrel{{\scriptstyle
BS2}}{{\rightarrow}}&\acute{t}\left|10\right\rangle_{D1,D2}+i\acute{r}\left|01\right\rangle_{D1,D2},\\\
\left|01\right\rangle_{bc}&\stackrel{{\scriptstyle
BS2}}{{\rightarrow}}&i\acute{r}\left|10\right\rangle_{D1,D2}+\acute{t}\left|01\right\rangle_{D1,D2},\end{array}$
(8)

where the real numbers $\acute{t}$ and $\acute{r}$ are the transmissivity and
reflectivity of the BS2, respectively, with $\acute{t}^{2}+\acute{r}^{2}=1$.
The output ports of the BS2 are $D1$ and $D2$. Putting Eq. (8) in Eq. (7), we
get the output state of BS2 as

$\left|Out\right\rangle_{BS2}=\frac{1}{\sqrt{2}}\left[\left(\acute{t}e^{i\theta}-i\acute{t}e^{i\theta}\varPhi\alpha
a^{\dagger}-\acute{r}\right)\left|\alpha\right\rangle_{a}\left|10\right\rangle_{D1,D2}+\left(i\acute{r}e^{i\theta}+i\acute{t}+\acute{r}e^{i\theta}\varPhi\alpha
a^{\dagger}\right)\left|\alpha\right\rangle_{a}\left|01\right\rangle_{D1,D2}\right].$
(9)

When the single photon is detected at $D1$, the coherent state
$\left|\alpha\right\rangle$ in mode $a$ is projected into the state

$\left|\varPsi\right\rangle=\frac{1}{\sqrt{2}}\left[\left(\acute{t}e^{i\theta}-\acute{r}\right)\left|\alpha\right\rangle_{a}-i\acute{t}e^{i\theta}\varPhi\alpha
a^{\dagger}\left|\alpha\right\rangle_{a}\right].$ (10)

We get the HCS, which is the linear superposition of the coherent state
$\left|\alpha\right\rangle$ and the single-photon-added coherent state
$a^{\dagger}\left|\alpha\right\rangle$ when $D1$ fires. The amplitude and
phase of the superposition of the CS and the SPAC state can be manipulated by
adjusting $\acute{r}$ or $\acute{t}$, $\varPhi$, $\alpha$ and $\theta$. In its
general form, we can write the HCS as

$\left|\varPsi\right\rangle=\mathscr{N}\left[\sqrt{\epsilon}\left|\alpha\right\rangle+\sqrt{1-\epsilon}e^{i\varphi}a^{\dagger}\left|\alpha\right\rangle\right],$
(11)

where $\epsilon$ determines the amplitudes of the CS and SPAC states in HCS
and $\varphi$ is the phase angle between the CS and the SPAC states. The
normalization constant $\mathcal{N}$ can be expressed as

$\mathscr{N}=\left[1+2\sqrt{\epsilon(1-\epsilon)}\,Re\left[\alpha
e^{-i\varphi}\right]+\left(1-\epsilon\right)\left|\alpha\right|^{2}\right]^{-\frac{1}{2}}.$
(12)

The parameter $\epsilon$ is real and has its value $0\leq\epsilon\leq 1$
determines the relative contribution of CS and the SPAC state in the HCS. For
$\epsilon=0$ and $\epsilon=1$, the HCS transforms to the SPAC state and the
CS, respectively. In what follows, the expression for the HCS in Eq. (11) is
exploited to investigate higher-order squeezing and higher-order antibunching.

##  4 Higher-order Squeezing

A nonclassical phenomenon known as squeezing involves reducing the uncertainty
in a particular quadrature at the cost of increasing the uncertainty in a
noncommutating quadrature. Higher-order squeezing is a notion in quantum
optics that extends the idea of squeezing past the uncertainty in the
quadrature measured via variance. Two separate approaches are typically used
to study higher-order squeezing [23, 24, 25]. Hong and Mandel originally
proposed an approach for defining higher-order squeezing in 1985 [23, 24],
which is a natural generalization of second-order squeezing. In contrast, the
second notion of higher-order squeezing was introduced in 1987 by Hillary
[25], where the reduction of the variance of the amplitude-powered quadrature
variable with respect to its coherent state counterpart was considered as
manifestation of higher-order squeezing. In this section, we will investigate
the higher-order Hong-Mandel-type squeezing of the hybrid coherent state. In
order to investigate the higher-order squeezing, we consider two quadrature
components of the field as $X_{\psi}$ and $X_{\psi+\frac{\pi}{2}}$. Let the
commutation relation between the quadratures is
$\left[X_{\psi},X_{\psi+\frac{\pi}{2}}\right]=iC$, where $C$ is a constant. In
the course of the work, we will assume that $\hbar=1.$ Now, the general
expression for the $2n$-order variance of any quadrature $X_{\psi}$ can be
written as [23]

$\begin{array}[]{lcl}\left\langle\left(\Delta
X_{\psi}\right)^{2n}\right\rangle&=&\left\langle\colon\left(\Delta
X_{\psi}\right)^{2n}\colon\right\rangle+\frac{2!}{1!}\binom{2n}{2}\frac{C}{4}\left\langle\colon\left(\Delta
X_{\psi}\right)^{2n-2}\colon\right\rangle\\\
&+&\frac{4!}{2!}\binom{2n}{4}\left(\frac{C}{4}\right)2\left\langle\colon\left(\Delta
X_{\psi}\right)^{2n-4}\colon\right\rangle+\cdots+\left(2n-1\right)!!\left(\frac{C}{2}\right)^{n}\\\
&=&S_{\psi}^{(2n)}+\left(2n-1\right)!!\left(\frac{C}{2}\right)^{n},\end{array}$
(13)

where the notation $::$ represents the normal ordering where the creation
operators are consistently maintained to the left of the annihilation
operators and

$S_{\psi}^{(2n)}=\sum_{m=0}^{n-1}\frac{(2n)!}{(2m+2)!(n-m-1)!}\left(\frac{C}{4}\right)^{n-m-1}\left\langle\colon\left(\Delta
X_{\psi}\right)^{2m+2}\colon\right\rangle.$ (14)

In the Eqs. (13) and (14), $m$ and $n$ are the positive integers. Now, the
falls of the $2n$-order fluctuation of a quadrature are below its CS value,
indicating the $2n$-order squeezing of that quadrature. So, the essential
requirement for $2n$-order Hong-Mandel type squeezing in $X_{\psi}$ quadrature
is [24, 26]

$\left\langle\left(\Delta
X_{\psi}\right)^{2n}\right\rangle<\left(2n-1\right)!!\left(\frac{C}{2}\right)^{n}.$
(15)

From the Eqs. (13), (14), and (15) it is clear that the signature of
$2n$-order squeezing of the quadrature $X_{\psi}$ will be seen if
$S_{\psi}^{(2n)}<0$. Now, for the present work, we define the quadrature
operators as
$X_{\psi}=\frac{1}{\sqrt{2}}\left(a^{\dagger}e^{i\psi}+ae^{-i\psi}\right)$ and
$X_{\psi+\frac{\pi}{2}}=\frac{i}{\sqrt{2}}\left(a^{\dagger}e^{i\psi}-ae^{-i\psi}\right)$.
Using the commutation relation defined above, we obtain the most general
expression of $S_{\psi}^{(2n)}$ as

$\begin{array}[]{lcl}S_{\psi}^{(2n)}&=&\sum_{m=0}^{n-1}\frac{(2n)!}{(n-m-1)!(2m+2)!}\left(\frac{C}{4}\right)^{(n-m-1)}\sum_{p=0}^{2m+2}\binom{2m+2}{p}\left(\frac{1}{\sqrt{2}}\right)^{2m+2}\\\
&\times&\sum_{l=0}^{2m+2-p}\binom{2m+2-p}{l}\left\langle
a^{\dagger(2m-p-l+2)}a^{l}\right\rangle
e^{i(2m-p-2l+2)\psi}\left(-1\right)^{p}\\\
&\times&\sum_{j=0}^{p}\binom{p}{j}\left\langle
a^{\dagger}\right\rangle^{p-j}\left\langle
a\right\rangle^{j}e^{i(p-2j)\psi}.\end{array}$ (16)

In our present study, $C=1.$ The expectation values of the $\left\langle
a^{\dagger n}a^{m}\right\rangle$ and $\left\langle a\right\rangle$ for our HCS
as defined in Eq. (11) are obtained in order to investigate the possibilities
of observing the higher-order squeezing. We derive

$\begin{array}[]{lcl}\left\langle a^{\dagger
n}a^{m}\right\rangle&=&\mathscr{N}^{2}\alpha^{*n-1}\alpha^{m-1}\left[\epsilon\left|\alpha\right|^{2}+\sqrt{\epsilon(1-\epsilon)}\left\\{e^{i\varphi}\left(m\alpha^{*}+\left|\alpha\right|^{2}\alpha^{*}\right)\right.\right.\\\
&&\left.\left.+e^{-i\varphi}\left(n\alpha+\left|\alpha\right|^{2}\alpha\right)\right\\}+\left(1-\epsilon\right)\left\\{nm+(n+m+1)\left|\alpha\right|^{2}+\left|\alpha\right|^{4}\right\\}\right],\end{array}$
(17)

$\left\langle
a\right\rangle=\mathscr{N}^{2}\left[\epsilon\alpha+\sqrt{\epsilon(1-\epsilon)}\left\\{e^{i\varphi}\left(1+\left|\alpha\right|^{2}\right)+e^{-i\varphi}\alpha^{2}\right\\}+\left(1-\epsilon\right)\alpha\left(2+\left|\alpha\right|^{2}\right)\right].$
(18)

Using Eqs. (17) and (18), we evaluate Eq. (16) to study the $2n$-order Hong-
Mandel type squeezing for HCS. In order to get the flavour of the
nonclassicality manifested through Eq. (16), we have plotted

(a)  Figure 2: (Colour online) Plot of higher-order squeezing with
$\left|\alpha\right|$: (a) in $X_{\psi}$quadrature for $n=1$, $2$, and $3$
with $\epsilon=0.5$; (b) in both quadratures for $\epsilon=0.5$.

the variation of $S_{\psi}^{(2n)}$ with $\left|\alpha\right|$ in Fig. 2. The
signature of higher-order squeezing with different orders is evident in Fig.
2(a), which also demonstrates that the degree of squeezing rises with order.
Fig. 2(b) shows that the squeezing in one quadrature prohibits the squeezing
in another quadrature as expected due to Heisenberg’s uncertainty relation.
Fig. 2 clearly shows that when $\left|\alpha\right|$ rises, the plot
approaches a coherent state. The solid line, dashed line, dotted line, and
dot-dashed line in Fig. 3

Figure 3: (colour online) Plot of $4^{th}$ order squeezing with
$\left|\alpha\right|$ for $\epsilon=0,\,0.25,\,0.5,$ and $0.75$.

represent the $4^{th}$-order squeezing for $\epsilon=0,$ $0.25,$ $0.5$, and
$0.75,$ respectively. Hence, the value of the mixing parameter $\epsilon$
determines the nature of the quadrature squeezing. As the
$\left|\alpha\right|$ increases, the HCS goes closer to the coherent state.
The HCS loses the nonclassical squeezing property for $\epsilon=1.$
Interestingly, for single photon ( $\left|\alpha\right|\rightarrow 0$) or few
photons ( $\left|\alpha\right|\rightarrow$small), the HCS ($\epsilon\neq 0$)
is observed to show higher-order squeezing, but is not observed if we take
only the SPAC state $(\epsilon=0)$. This is the advantageof considering the
HCS over the SPAC state.

##  5 Higher-order antibunching

Photon antibunching is another well-known nonclassical effect of our interest.
Photon antibunching is a phenomenon where photons prefer to be discrete rather
than clustered. In this section, we investigate the higher-order as well as
the normal-order antibunching properties of the HCS. Lee [27] first introduced
specific criteria for higher-order antibunching. Numerous, nearly identical
variations of this criterion have been offered over time. Pathak and Garcia
presented one such criterion [28]. Higher-order antibunching in various
optical systems [29, 30] and atomic systems [31, 32] is extensively studied by
some of the present authors using the Pathak and Garcia criterion. But no
effort has been made yet to investigate the existence of higher-order
antibunching in HCS. In order to do this, we use the Pathak-Garcia criterion
[28],which is for $n^{th}$order antibunching as follows:

$g^{(n+1)}=\frac{\left\langle
a^{\dagger\,n+1}a^{n+1}\right\rangle}{\left\langle
a^{\dagger}a\right\rangle^{n+1}}<1,$ (19)

where $n=1$ and $n\geq$2 refer to usual antibunching, and higher-order
antibunching, respectively. Now, in order to investigate the higher-order
antibunching we obtain $\left\langle a^{\dagger}a\right\rangle$ for the HCS as

$\left\langle
a^{\dagger}a\right\rangle=\mathscr{N}^{2}\left[1-\epsilon+\left(3-2\epsilon\right)\left|\alpha\right|^{2}+\left(1-\epsilon\right)\left|\alpha\right|^{4}+2\sqrt{\epsilon(1-\epsilon)}\left(1+\left|\alpha\right|^{2}\right)\,Re\left[\alpha
e^{-i\varphi}\right]\right].$ (20)

Using Eqs. (17) and (20) we evaluate Eq. (19) to study the higher-order
antibunching of the HCS. To gain a sense of the higher-order antibunching in
HCS, we plot Eq. (19) with $\left|\alpha\right|$ in Fig. 4 which shows that
antibunching is observed in normal as well as higher order. Fig. 4(a) shows
that the degree of nonclassicality rises with order number. Fig. 4(b) shows
how the antibunching property also depends on the mixing parameter $\epsilon$
and it reduces into a coherent state distribution for $\epsilon=1$.

Figure 4: (Colour online) Plot of higher-order antibunching with
$\left|\alpha\right|$: (a) for $n=1,$$2$ and $3$ with $\epsilon=0.5$; (b) for
$\epsilon=0,\,0.5,\,1,$with $n=2$.

##  6 Conclusions

Nonclassical properties and/or quantum states are needed for obtaining any
kind of quantum advantage. In fact, it’s essential for achieving quantum
supremacy [33]. Further, an important resource for achieving useful
nonclassical features is the set of engineered quantum states. Several
applications of engineered quantum states have been proposed in the recent
past, but not all such states can be produced using the available
technologies. This makes the experimentally realizable engineered quantum
states having useful nonclassicality extremely important. The importance of
the present work is here, as the present study shows that HCS is an
experimentally realizable engineered quantum state with higher-order
nonclassical features. Although up to the $6^{th}$order squeezing and
$4^{th}$order antibunching have been illustrated, our results hold for any
order. It’s interesting to note that the HCS exhibits higher-order squeezing
for low photon counts, whereas the SPAC state does not. This illustrates why
researching the HCS is more important than the SPAC state. Additionally, the
depth of antibunching increases with increasing orders, indicating the
necessity of studying higher-order antibunching.

As the very clear presence of higher-order antibunching and higher-order
squeezing in HCS is obtained here following an analytical approach, it’s
expected that the beauty of the analytic solution will encourage others to
perform similar studies with other engineered quantum states. Further, HCS,
being higher-order nonclassical and experimentally realizable, is expected to
be useful in quantum metrology and quantum communication in the near future.
With this hope, we conclude this paper.

Acknowledgement: B.S. expresses gratitude to Anirban Pathak for some valuable
technical conversations with him and his assistance and interest in this work.

## References

  * [1] A. Zavatta, S. Viciani, and M. Bellini, Phys. Rev. A 72 (2005) 023820.
  * [2] E. Bimbard, N. Jain, A. MacRae, and A. I. Lvovsky, Nat. Photonics 4 (2010) 243.
  * [3] P. M. Harrington, E. J. Mueller, and K. W. Murch, Nat. Rev. Phys. 4 (2022) 660.
  * [4] G. S. Agarwal and K. Tara, Phys. Rev. A 43 (1991) 492.
  * [5] R. Román-Ancheyta, C. G. Gutiérrez, and J. Récamier, J. Opt. Soc. Am. B 31 (2014) 38.
  * [6] A. Zavatta, S. Viciani, and M. Bellini, Science 306 (2004) 660 .
  * [7] S. Sivakumar, Phys. Rev. A 83 (2011) 035802.
  * [8] M. Miranda and D. Mundarain, Quantum Inf. Process. 16 (2017) 298.
  * [9] D. Wang, M. Li, F. Zhu, Z-Q. Yin, W. Chen, Z-F. Han, G-C. Guo, and Q.Wang, Phys. Rev. A 90 (2014) 062315.
  * [10] P. V. P. Pinheiro and R. V. Ramos, Quantum Inf. Process. 12 (2013) 537.
  * [11] J. T. Francis and M. S. Tame, Phys. Rev. A 102 (2020) 043709.
  * [12] Y. Turek, A. Islam, and A. Abliz, Eur. Phys. J. Plus 138 (2023)72.
  * [13] W. J. Xu, T. Yusufu, and Y. Turek, Phys. Rev. A 105 (2022) 022210.
  * [14] A. I. Lvovsky, H. Hansen, T. Aichele, O. Benson, J. Mlynek, and S. Schiller, Phys. Rev. Lett. 87 (2001) 050402.
  * [15] Y. Turek, N. Aishan, and A. Islam, Phys. Scr. 98 (2023) 075103.
  * [16] T. M. Duc and J. Noh, Opt. Commun. 281 (2008) 2842.
  * [17] Z. Li, H. Guo, H. Liu, J. Li, H. Sun, R. Yang, K. Liu, and J. Gao, Adv. Quantum Technol. 5 (2022) 2200055.
  * [18] S. Pirandola, U. L. Andersen, L. Banchi, M. Berta, D. Bunandar, R. Colbeck, D. Englund, T. Gehring, C. Lupo, C. Ottaviani, J. L. Pereira, M. Razavi, J. S. Shaari, M. Tomamichel, V. C. Usenko, G. Vallone, P. Villoresi, and P. Wallden, Adv. Opt. Photonics 12 (2020) 1012.
  * [19] K. Laiho, T. Dirmeier, M. Schmidt, S. Reitzenstein, and C. Marquardt, Phys. Lett. A, 435 (2022) 128059\.
  * [20] A. Pathak and A. Ghatak, J. Electromagn. Waves Appl. 32 (2018) 229.
  * [21] B. M. Escher, A. T. Avelar, T. M. da Rocha Filho, and B. Baseia, Phys. Rev. A. 70 (2004) 025801.
  * [22] F. Matsuoka, A. Tomita, and Y. Shikano, Quantum Sutd.: Math. Found. 4 (2017) 159.
  * [23] C. K. Hong and L. Mandel, Phys. Rev. Lett. 54 (1985) 323.
  * [24] C. K. Hong and L. Mandel, Phys. Rev. A 32 (1985) 974.
  * [25] M. Hillery, Phys. Rev. A 36 (1987) 3796.
  * [26] A. Verma and A. Pathak, (2010) arXiv:1004.1689v1 [quant-ph].
  * [27] C. T. Lee, Phys. Rev. A 41 (1990) 1721\.
  * [28] A. Pathak, and M.E. Garcia, Appl. Phys. B 84 (2006) 479.
  * [29] K. Thapliyal, A. Pathak, B. Sen, and J.Peřina, Opt. Commun. 444 (2019) 111.
  * [30] K. Thapliyal, A. Pathak, B. Sen, and J. Peřina, Phys. Rev. A 90 (2014) 013808.
  * [31] S. K. Giri, B. Sen, C. H. R. Ooi, and A. Pathak, Phys. Rev. A 89 (2014) 033628.
  * [32] S. K. Giri, K. Thapliyal, B. Sen, and A. Pathak, Physica A 466 (2017) 140.
  * [33] T. Douce, D. Markham, E. Kashefi, E. Diamanti, T. Coudreau, P. Milman, P. Van Loock, and G. Ferrini, Phys. Rev. Lett. 118 (2017) 070503.


# Strongly disordered Anderson insulator chains with generic two-body
interaction

B. Krajewski  Institute of Theoretical Physics, Faculty of Fundamental
Problems of Technology,
Wrocław University of Science and Technology, 50-370 Wrocław, Poland  L.
Vidmar  Department of Theoretical Physics, J. Stefan Institute, SI-1000
Ljubljana, Slovenia  Department of Physics, Faculty of Mathematics and
Physics, University of Ljubljana, SI-1000 Ljubljana, Slovenia  J. Bonča
Department of Physics, Faculty of Mathematics and Physics, University of
Ljubljana, SI-1000 Ljubljana, Slovenia  Department of Theoretical Physics, J.
Stefan Institute, SI-1000 Ljubljana, Slovenia  M. Mierzejewski  Institute of
Theoretical Physics, Faculty of Fundamental Problems of Technology,
Wrocław University of Science and Technology, 50-370 Wrocław, Poland

(August 28, 2024)

###### Abstract

The random-field spin-1/2 XXZ chains, and the corresponding Anderson
insulators of spinless fermions with density-density interaction, have been
intensively studied in the context of many-body localization. However, we
recently argued [Phys. Rev. Lett. 129, 260601 (2022)] that the two-body
density-density interaction in these models is not generic since only a small
fraction of this interaction represents a true local perturbation to the
Anderson insulator. Here we study ergodicity of strongly disordered Anderson
insulator chains choosing other forms of the two-body interaction for which
the strength of the true perturbation is of the same order of magnitude as the
bare two-body interaction. Focusing on the strong interaction regime,
numerical results for the level statistics and the eigenstate thermalization
hypothesis are consistent with emergence of ergodicity at arbitrary strong
disorder.

##  I Introduction

Extensive numerical studies concerning ergodicity of disorders systems have
been carried out mostly for the one-dimensional random-field spin-1/2 XXZ
model (see, e.g., pal10 ; barisic10 ; luitz15 ; luitz16 ; torres15 ; bera15 ;
Hauschild_2016 ; Devakul2015 ; bertrand_garcia_16 ; Doggen2018 ;
suntajs_bonca_20a ; suntajs_2020 ; sierant_lewenstein_20 ). This spin model
can be mapped onto a chain of spinless fermions, for which the Hamiltonian can
be expressed as

$H=H_{0}(W)+H^{\prime}(V)\;,$ (1)

where the single-particle term $H_{0}(W)$ describes the Anderson insulator
with the disorder strength $W$ and the hopping amplitude set to unity, while
the interaction $H^{\prime}(V)$ is given by the density-density interaction of
strength $V$. This is one of the simplest models which allows to study the
influence of two-body interaction on the Anderson localization.

At a sufficiently strong disorder, $W>W^{*}$, numerical studies of finite
systems indicated that localization may also persist in the presence of the
two-body interaction, leading to the concept of many-body localization (MBL)
basko06 ; oganesyan07 . In particular, the presence of MBL has been suggested
based on extremely slow dynamics znidaric08 ; serbyn13_1 ; mierzejewski2016 ;
luschen17 ; bordia2017_1 ; serbyn2017 ; Bera2019 ; chanda2020 , suppressed
transport berkelbach10 ; lev15 ; barisic16 ; steinigeweg16 ; bera2017 ;
prelovsek116 , entanglement entropy properties znidaric08 ; bardarson12 ;
kjall14 ; baygan15 ; pietracaprina_parisi_17  and level statistics analyses
oganesyan07 ; pal10 ; luitz15 ; zakrzewski16 ; sierant_zakrzewski_19 .

In spite of these results, however, the fate of localization in macroscopic
systems remains an open problem. In particular, the transition/crossover
disorder $W^{*}$, as estimated from the energy spectrum suntajs_bonca_20a ;
suntajs_2020 , shows a linear drift with the system size $L$, namely
$W^{*}\propto L+{\rm const}$. The latter observations as well as several other
results obtained from subsequent studies suggest that macroscopic systems may
be ergodic kieferemmanouilidis_unanyan_20 ; sels2020 ;
kieferemmanouilidis_unanyan_21 ; leblond_sels_21 ; vidmar2021 ;
Sels_dilute_2021 ; Sels_2022 ; evers2023  and diffusive, although with a very
small diffusion constant barisic10 . At even larger disorders, much larger
than the previous estimates of $W^{*}$ sierant2022 ; Morningstar2022 , the
state of the art numerical approaches may not provide a conclusive answer to
this questionabanin2019 ; Panda2020 ; abanin_bardarson_21 ;
sierant_lewenstein_20 ; crowley_chandran_22 ; bolther_kehrein_22 .

A nontrivial feature of the $V-W$ state diagram in finite systems is the
emergence of reentrant behavior deluca13 ; lev15 . This means that at fixed
$W$, ergodicity is restored at some nonzero interaction $V_{1}$, while at some
larger interaction $V_{2}>V_{1}$ one again observes the breakdown of
ergodicity. This property may suggest that the density-density interaction
cannot be considered as a perturbation that just destroys the Anderson
localization of noninteracting particles but it may also stabilize
localization once it is large enough.

In a recent Letter krajewski_vidmar_22 , we showed that only a small part of
the density-density interaction $H^{\prime}(V)$, termed true local
perturbation (shortly, true perturbation), does not commute with $H_{0}$.
Namely, the major part of $H^{\prime}(V)$ can be expressed via products of
occupations of the single-particle Anderson states, and hence this part of
$H^{\prime}(V)$ does not perturb the Anderson insulator. At strong disorder,
the true perturbation may hence be too small to be correctly captured by
numerical calculations carried out in finite systems. By rescaling the true
perturbation so that its strength is of the same order as the other terms of
the Hamiltonian, we introduced a rescaled model that exhibits ergodicity at
arbitrary strong $W$ krajewski_vidmar_22 . However, this approach is
computationally demanding and the spatial structure of the true perturbation
is rather complex.

In this work, we demonstrate that the smallness of the true perturbation is
specific for models with the density-density interaction, and that other forms
of the two-body interaction do not share this property. Therefore, the
emergence and breakdown of ergodicity in the latter systems can be studied
without applying a complex procedure that singles out and rescales the true
perturbation. As the main result, we show evidence that a sufficiently strong
two-body interaction restores ergodicity even at extremely large disorder such
as $W\sim 20$. This observation also suggests the absence of reentrant
behavior at fixed disorder $W$ when interaction is increased.

The paper is organized as follows. In Sec. II, we recall the notion of the
true perturbation and introduce the studied forms of the two-body
interactions. Then, in Sec. III we discuss numerical results obtained for the
level statistics and the eigenstate thermalization hypothesis (ETH) deutsch_91
; srednicki_94 ; rigol_dunjko_08 ; dalessio_kafri_16 . We summarize and
discuss our results in Sec. IV.

##  II Models and the strength of the true perturbation

We consider a chain of length $L$ containing $L/2$ spinless fermions, where we
assume periodic boundary conditions. We discuss the ergodicity of three
distinct models of the form (1), which share the same single-particle part,
$H_{0}$, but have different two-body interaction terms, $H^{\prime}$. In
particular, we take

$H_{0}=\frac{t_{h}}{2}\sum_{i}\left(c^{\dagger}_{i+1}c_{i}+{\rm
H.c.}\right)+\sum_{i}\epsilon_{i}\left(n_{i}-\frac{1}{2}\right),$ (2)

where $c^{\dagger}_{i}$ creates a fermion at site $i$,
$n_{i}=c^{\dagger}_{i}c_{i}$ and $\epsilon_{i}$ is a random potential (i.e.,
the disorder) uniformly distributed in the box $\epsilon_{i}\in[-W,W]$. From
this point on, we set $t_{h}=1$. We refer to the model in Eq. (2) the Anderson
model. The noninteracting part is diagonal in the Anderson single-particle
basis:

$H_{0}=\sum_{\alpha}\frac{1}{2}\varepsilon_{\alpha}Q_{\alpha}+{\rm
const},\quad\quad Q_{\alpha}=2n_{\alpha}-1,$ (3)

where $n_{\alpha}$ is the particle number operator for the Anderson state
$|\alpha\rangle$ with the energy $\varepsilon_{\alpha}$, and $Q_{\alpha}$ are
referred to as the Anderson local integrals of motions (shortly, Anderson
LIOMs). It is convenient to single out the operator density, $h(i)$, of the
two-body interaction, and express $H^{\prime}$ as

$H^{\prime}=V\sum_{i}h(i)\;.$ (4)

The first model we have in mind is the standard model of MBL with the nearest-
neighbor density-density (dd) interaction for spinless fermions,

$h_{\rm
dd}(i)=\left(n_{i}-\frac{1}{2}\right)\left(n_{i+1}-\frac{1}{2}\right)\;,$ (5)

which maps onto the $S^{z}_{i}S^{z}_{i+1}$ interaction in the random-field
spin-1/2 XXZ model. Increasing the disorder in the Anderson model in Eq. (2)
yields a shorter localization length, and ultimately at large $W$ each single-
particle Anderson wave function $\langle i|\alpha\rangle$ is strongly peaked
at a single site $i=i(\alpha)$. Then, the occupations in Eq. (5) become
similar to the Anderson LIOMs in Eq. (3),
$n_{i(\alpha)}-\frac{1}{2}\simeq\frac{1}{2}Q_{\alpha}$. Therefore, the value
of $V$ alone does not specify whether $H^{\prime}_{\rm dd}=V\sum_{i}h_{\rm
dd}(i)$ represents strong or weak perturbation to $H_{0}$.

In order to quantify this intuitive understanding [especially the similarity
between $n_{i(\alpha)}-\frac{1}{2}$ and $Q_{\alpha}$], we use the Hilbert-
Schmidt inner product of two operators $\langle AB\rangle=\frac{1}{Z}{\rm
Tr}(AB)$ and the corresponding norm of operators $||A||^{2}=\langle
AA\rangle$. Here, the trace is carried out over the many-body states and
$Z={\binom{L}{L/2}}$ is the dimension of the many-body Hilbert space. We split
the perturbation operators from Eq. (5) into two parts,

$h_{\rm dd}(i)=h^{\parallel}_{\rm dd}(i)+h^{\perp}_{\rm dd}(i)\;,$ (6)

in such a way that both parts are mutually orthogonal $\langle
h^{\parallel}_{\rm dd}(i)h^{\perp}_{\rm dd}(i)\rangle=0$ and
$h^{\parallel}_{\rm dd}(i)$ commutes with $H_{0}$. Therefore, only the latter
part, $Vh^{\perp}_{\rm dd}(i)$, represents the density of the true local
perturbation, and the operators can be expressed as (see krajewski_vidmar_22
for details)

$\displaystyle h^{\parallel}_{\rm dd}(i)$ $\displaystyle=$
$\displaystyle\sum_{\alpha,d\leq d_{\rm max}}\langle h_{\rm
dd}(i)Q_{\alpha}Q_{\alpha+d}\rangle\;\;Q_{\alpha}Q_{\alpha+d}\;,$ (7)
$\displaystyle h^{\perp}_{\rm dd}(i)$ $\displaystyle=$ $\displaystyle h_{\rm
dd}(i)-h^{\parallel}_{\rm dd}(i)\;.$ (8)

The index $\alpha+d$ of $Q_{\alpha+d}$ in Eq. (7) corresponds to the single-
particle Anderson state for which the maximum of the wave function $|\langle
i|\alpha+d\rangle|$ is shifted left by $d$ sites with respect to the maximum
of $|\langle i|\alpha\rangle|$. We recall that at strong disorder the strength
of the true perturbation decays with $W$ as $V||h^{\perp}_{\rm dd}(i)||\sim
V/W$, see Fig. 1 in krajewski_vidmar_22 , whereas $Vh^{\parallel}_{\rm
dd}(i)\sim V$ should be considered as a part of the unperturbed energy
density.

As a consequence of the above analysis, the strongly disordered XXZ models
always represent weakly perturbed Anderson insulators, independently of $V$.
Weak perturbations are challenging for numerical calculations which are
carried out in finite systems. To avoid possible numerical artifacts, in Ref.
krajewski_vidmar_22  we rescaled the strength of the true perturbation in such
a way that its norm is of the same order of magnitude as the norm of $H_{0}$.
Although the rescaled model was shown to be ergodic, the explicit form of the
rescaled Hamitonian is rather complex and may raise concerns whether the
conclusions concerning its ergodicity are generic.

To avoid difficulties originating from the presence of a very weak
perturbation, in this work we study simple Hamiltonians in which the norms of
the single-particle term and the true two-body perturbation are controlled by
two independent parameters $W$ and $V$, respectively. Since the smallness of
the true perturbation in the XXZ chains is specific for the density-density
interaction, in this work, we study systems with two-body interactions that
cannot be expressed in terms of occupations of the lattice sites. We mostly
focus on the two-body interaction $h(i)$ from Eq. (4), which is given by a
pair-hopping (ph) term

$h_{\rm ph}(i)=c^{\dagger}_{i-1}c^{\dagger}_{i}c_{i+1}c_{i+2}+{\rm
H.c.}\;\;\;.$ (9)

In order to demonstrate that our conclusions do not originate from a specific
choice of the perturbation, at the end of the paper we also discuss numerical
results obtained for the operator

$\tilde{h}_{\rm ph}(i)=c^{\dagger}_{i}c^{\dagger}_{i+1}c_{i+2}c_{i-1}+{\rm
H.c.}\;\;\;.$ (10)

Results obtained for both models are qualitatively the same. Namely, we show
that even at a very strong disorder $W$, a sufficiently strong perturbation
restores the ergodicity of the studied model. The larger the system, the
weaker the perturbation for which the ergodicity is restored.

##  III Numerical results

Figure 1: Norms of the operators $||h_{\rm ph}^{\parallel}(i)||$ defined via
Eq. (7) by replacing $h_{\rm dd}(i)\to h_{\rm ph}(i)$. (a) Probability density
functions (PDFs) of $||h_{\rm ph}^{\parallel}(i)||$ at $L=14$ and $d_{\rm
max}=2$ for various $W$. The distributions are calculated at all sites for
2000 realizations of the disorder. (b) $W$-dependence of the medians of the
distributions shown in (a). Dashed line is the function $\propto 1/W^{3.5}$.

###  III.1 The strength of the true perturbation

We now focus on the Anderson insulators with pair-hopping interaction,

$H=H_{0}+H^{\prime}(V)=H_{0}+V\sum_{i}h_{\rm ph}(i),$ (11)

where $h_{\rm ph}(i)$ is defined in Eq. (9). First, we demonstrate that the
strength of the true perturbation in this model only weakly depends on $W$ and
is of the same order as the strength of the bare two-body interaction,
$V||h_{\rm ph}(i)||$. To this end, we replace $h_{\rm dd}(i)\to h_{\rm ph}(i)$
in Eq. (6) and we use the projections defined via Eqs. (7) and (8) to
calculate $h^{\parallel}_{\rm ph}(i)$ and $h^{\perp}_{\rm ph}(i)$,
respectively. Due to the orthogonality of these two terms one obtains the
following identity:

$||h_{\rm ph}(i)||^{2}=||h^{\parallel}_{\rm ph}(i)||^{2}+||h^{\perp}_{\rm
ph}(i)||^{2},$ (12)

where the norm on the left-hand side of Eq. (12) does not depend on $W$ or
$V$. Therefore, the strength of the true perturbation, $V||h^{\perp}_{\rm
ph}(i)||$, can be estimated indirectly from the component which commutes with
$H_{0}$, i.e., from $h^{\parallel}_{\rm ph}(i)$.

Using exact diagonalization for a system with $L=14$ sites, we calculated the
norms of $||h^{\parallel}_{\rm ph}(i)||$ for various realizations of disorder
and various sites $i$. From the latter we determined the probability density
functions (PDFs) for $||h^{\parallel}_{\rm ph}(i)||$, which are shown in Fig.
1(a) for various $W$. In Fig. 1(b) we show medians of these distributions as a
function of the disorder strength $W$. Upon increasing $W$, the maxima of
these distributions shift towards smaller values of $||h^{\parallel}_{\rm
ph}(i)||$. Moreover, the medians decrease faster than $1/W^{3}$. This effect
is opposite to the results for the density-density interactions, where
$||h^{\parallel}_{\rm dd}(i)||$ are very close to $|||h_{\rm dd}(i)||$,
whereas $|||h^{\perp}_{\rm dd}(i)||$ decrease as $1/W$. In the present case
$||h^{\parallel}_{\rm ph}(i)||$ are negligible, in particular at strong
disorder.

To summarize, based on the results in Fig. 1 and using the identity from Eq.
$(\ref{norms})$, we conclude that the strength of the true perturbation in the
model (11) is determined mainly by $V$ and, up to reasonable accuracy, one may
assume that it is equal to the bare interaction $H^{\prime}$.

###  III.2 Average level spacing ratio

To identify the ergodic and nonergodic regimes of the model from Eq. (11), we
start with the commonly studied ergodicity indicator, i.e., the ratio
$r_{n}={\rm min}(\delta E_{n},\delta E_{n+1})/{\rm min}(\delta E_{n},\delta
E_{n+1})$ of the nearest level spacings $\delta E_{n}=E_{n+1}-E_{n}$, where
$E_{n}$ are (sorted) energy levels oganesyan07 ; atas2013 . We average $r_{n}$
over the middle third of the energy spectrum, as well as over $4000$
realizations of the disorder, and denote the average level spacing ratio as
$\langle r\rangle$.

Figures 2(a), 2(c) and 2(e) show results for $\langle r\rangle$ versus $L$ at
$V=1,3,6$, respectively, obtained at $W=2,3,4,5,6,7,8,10$. In contrast, Figs.
2(b), 2(d) and 2(f) show results for $\langle r\rangle$ versus $W$ and the
same $V=1,3,6$, respectively, obtained at $L=10,12,14,16$.

Figure 2: Average level spacing ratio $\langle r\rangle$ for the model in Eq.
(11). The averaging is carried out over $Z/3$ levels from the middle of the
spectrum and over 4000 realizations of disorder. (a,c,e) Dependence of
$\langle r\rangle$ on $L$ at interactions $V=1,3,6$, respectively. Dark to
light colors denote the disorders $W=2,3,4,5,6,7,8,10$, see also the colorbar
in (e). (b,d,f) Dependence of $\langle r\rangle$ on $W$ at interactions
$V=1,3,6$, respectively. Dark to light colors denote the system sizes
$L=10,12,14,16$, see also the legend in (f).  Figure 3: Average level spacing
ratio $\langle r\rangle$ for the model in Eq. (11). (a,b) Heat maps of
$\langle r\rangle$ at $L=12$ and 16, respectively. The solid line denotes
$V=W$. (c) Dependence of $\langle r\rangle$ on $L$ at $V=W$. Dark to light
colors denote the disorders $W=1,2,3,4,5,6,7,8,10,12,14,16,20,22,24$, see also
the colorbar. (d) Dependence of $\langle r\rangle$ on $W/L$ (main panel) and
$W/L^{2}$ (inset) at $V=W$. Dark to light colors denote the system sizes
$L=10,12,14,16,18$, see also the legend.

The results at weak perturbation $V=1$, see Figs. 2(a) and 2(b), share some
similarities with the numerical results obtained previously for disordered
chains with density-density interaction luitz15 ; bertrand_garcia_16 ;
Husex2017 ; suntajs_2020 ; sierant_prl_2020 . In particular, at sufficiently
large disorder, see Fig. 2(a), $\langle r\rangle$ decreases with $L$ towards
the value for Poisson level statistics $r_{\rm Pois}\simeq 0.39$, eventually
indicating nonergodicity. The finite-size phenomenology of the results for
$\langle r\rangle$ in systems with density-density interaction were
extensively studied in the past, giving rise to conjectures either about
robustness of ergodicity suntajs_2020  or about the convergence to a critical
point sierant_prl_2020 . Our study does not contribute significantly to new
results in this regard.

In contrast to systems with density-density interaction, in the present model
(11) one may tune the strength of the interaction and explore a broader regime
of parameters where the numerical results do not suffer from suppression of
the true perturbation. Figures 2(c)-2(f) show results for $\langle r\rangle$
at stronger interactions $V=3$ and 6. They show a strong tendency towards
restoring ergodicity. For example, in the system with $L=16$ sites, it is
reasonable to expect $\langle r\rangle$ to approach Gaussian orthogonal
ensemble (GOE) prediction $r_{\rm GOE}\simeq 0.53$ at disorders of strengths
$W\approx 10$. Moreover, the flow of the results with increasing $L$ suggests
restoration of ergodicity also at disorders much larger than $W\approx 10$.

A different perspective of the finite-size phenomenology of $\langle r\rangle$
is presented in Fig. 3. The values of $\langle r\rangle$ in the $V-W$ plane of
parameters are shown as heat maps in Figs. 3(a) and 3(b) at $L=12$ and 16,
respectively. In both cases, signatures of a reentrant behavior may be
observed, i.e., at a fixed $W$ and by increasing $V$, $\langle r\rangle$ first
increases towards $r_{\rm GOE}$ and then decreases again. However, another
important observation is also that by increasing the system size from $L=12$
to $L=16$, there is a strong enhancement of the region in which $r\approx
r_{\rm GOE}$.

To study the fate of the reentrant behavior when increasing the system size,
we focus on the parameter line $V=W$, which is shown as a solid line in Figs.
3(a) and 3(b). If the system exhibits a tendency towards restoring ergodicity
along this line, it is unlikely that the reentrant behavior survives in the
thermodynamic limit.

Figure 3(c) shows results for $\langle r\rangle$ versus $L$ at $V=W$. Even for
the largest disorder under consideration, $W=24$, we clearly observe a drift
of $\langle r\rangle$ towards $r_{\rm GOE}$ upon increasing $L$. A different
view on these results can be obtained from Fig. 3(d), in which we show
$\langle r\rangle$ versus $W/L$ at $V=W$. If the characteristic disorder at
which $\langle r\rangle$ starts to depart from $r_{\rm GOE}$ drifts linearly
with $L$, as in the case of density-density interactions suntajs_2020 , one
should observe a collapse of curves for $\langle r\rangle$ at different $L$
when plotted versus $W/L$. However, the results in Fig. 3(d) show signatures
that the drift is faster than linear. In particular, the values of disorder at
which $\langle r\rangle$ starts to deviate from $r_{\rm GOE}$ seem to increase
with system size as $W^{*}\propto L^{2}$, see the inset in Fig. 3(d). The
latter scaling is different from the scaling $W^{*}\propto L$ at $V=W$ line
for the density-density interaction (not shown). We interpret this result as a
strong evidence that the system in Eq. (11) becomes ergodic along the $V=W$
parameter line in the thermodynamic limit; hence the state diagram is not
expected to exhibit reentrant behavior.

###  III.3 Diagonal matrix elements of observables and the ETH

Next, we demonstrate the tendency of the studied model (11) to become ergodic
via testing the ETH. We perform numerical calculations for the site occupation
operator $A_{i}=2n_{i}-1$. This observable (or some linear combinations of
$A_{i}$) was commonly studied in the context of MBL mierzejewski2016 ; luitz16
; sirker14 ; pal10  because localization would imply that the matrix elements
of $A_{i}$ must violate the ETH. The latter violation would show up in the
nonvanishing fluctuations of the diagonal matrix elements $\langle
E_{m}|A_{i}|E_{m}\rangle$, which determine the value of the infinite-time
correlation function $\lim_{t\to\infty}\langle A_{i}(t)A_{i}\rangle$.

In order to estimate the fluctuations of the diagonal matrix elements, we
calculate the average eigenstate-to-eigenstate fluctuations Kim_strong2014 ;
Mondaini2016 ; jansen_stolpp_19 ,

$\langle\delta A_{i}\rangle=\frac{1}{\cal{Z}}\sum_{m}|\langle
E_{m+1}|A_{i}|E_{m+1}\rangle-\langle E_{m}|A_{i}|E_{m}\rangle|\;,$ (13)

where the eigenstate-to-eigenstate fluctuations are averaged over
${\cal{Z}}=Z/5$ eigenstates from the middle of the many-body spectrum.

Figure 4: Fluctuations of the diagonal matrix elements for the model defined
in Eq. (11). (a) and (b): PDFs of $\langle\delta A_{i}\rangle$, see Eq. (13),
at $V=6$ for a single site $i$ and 4000 realizations of disorder. (c)-(f):
Medians of distributions as in (a) and (b) versus $L$ for various $V$ and $W$.

Figures 4(a) and 4(b) show the PDFs of $\langle\delta A_{i}\rangle$ obtained
from a collection of results for 4000 realizations of disorder. In the case of
a finite system, $A_{i}$ at different sites are not independent since
$\sum_{i}A_{i}=0$. Therefore, we calculate $A_{i}$ for a single lattice site
$i$. These PDFs were obtained at $V=6$, when the average level spacing ratio
$\langle r\rangle$ matches the GOE prediction. At strong disorder $W=10$, see
Fig. 4(b), the values of $\langle\delta A_{i}\rangle$ remain large and the
presented distributions are broad. The latter feature indicates the presence
of substantial sample-to-sample fluctuation, which seems to be unavoidable in
small systems at strong disorder mierzejewski2020 ; krajewski2022 .
Nevertheless, upon increasing $L$ the distributions become narrower and their
maxima shift towards smaller eigenstate-to-eigenstate fluctuations. This
indicates that the diagonal matrix elements of $A_{i}$ satisfy the ETH.

In order to study in more detail how the distributions of $\langle\delta
A_{i}\rangle$ depend on $V$ and $W$, in Figs. 4(c)-4(f) we show the medians of
these distributions versus $L$. The conclusions are quite similar to those
derived previously from the level statistics. At weak two-body interactions,
see the results for $V=1$ in Figs. 4(e) and 4(f), strongly disordered systems
seem to violate the ETH since the eigenstate-to-eigenstate fluctuations stay
large for all accessible system sizes. However, the $L$-dependence of the
medians at $V>1$ exhibits clear negative curvatures in this regime of
parameters, suggesting that the eigenstate-to-eigenstate fluctuations should
eventually start decaying for sufficiently large systems. In the case of
stronger two-body interaction, which is the main focus of our study, such
decay is clearly visible for all studied disorder strengths.

###  III.4 Results for two-body interaction from Eq. (10)

Figure 5: Ergodicity indicators for the model defined in Eq. (14). (a) and
(b): Average level spacing ratio $\langle r\rangle$ versus $L$ at $V=3$ and 6,
respectively. Results are analogous to those in Figs. 2(c) and 2(e),
respectively, shown for the model from Eq. (11). (c) and (d) Medians of the
PDFs of the eigenstate-to-eigenstate fluctuations $\langle\delta A_{i}\rangle$
from Eq. (13) at $W=2$ and 10, respectively. Results are analogous to those in
Figs. 4(c) and 4(f), respectively, shown for the model from Eq. (11).

To demonstrate that the ergodicity of the model studied so far is not
restricted to the particular form of the interaction, in this section we
discuss the average level spacing ratio and the fluctuations of the diagonal
matrix elements of observables for the Hamiltonian

$\tilde{H}=H_{0}+V\sum_{i}\tilde{h}_{\rm ph}(i),$ (14)

where the pair-hoping term, $\tilde{h}_{\rm ph}(i)$, is given by Eq. (10).
Similarly to the previously discussed case also $||\tilde{h}^{\parallel}_{\rm
ph}(i)||$ decays with increasing $W$ (not shown); therefore, the entire two-
body interaction can be considered as a perturbation to $H_{0}$.

The average level spacing ratios $\langle r\rangle$ obtained at weak
interaction $V=1$ (not shown) are very similar to the results for the
previously studied model from Eq. (11) shown in Figs. 2(a) and 2(b). The
results at stronger interactions $V=3$ and 6 are shown in Figs. 5(a) and 5(b),
respectively. Although deviations of $\langle r\rangle$ from $r_{\rm GOE}$ are
more pronounced than in the case of the model from Eq. (11), the
$L$-dependence of $\langle r\rangle$ consistently supports that the spectral
properties of macroscopic systems should approach the GOE predictions.

Finally, we also calculated the fluctuations of the diagonal matrix elements
$\langle\delta A_{i}\rangle$ from Eq. (13). Similarly to the discussion in the
Sec. III.3, we obtained the PDFs of $\langle\delta A_{i}\rangle$, and in Figs.
5(c) and 5(d) we show the corresponding medians. Accessible system sizes do
not allow one to formulate univocal claims concerning ergodicity of
macroscopic systems at very strong disorder and weak two-body interactions.
However, results in the strong interaction regime clearly indicate that the
system satisfies the ETH.

##  IV Discussion

We studied the ergodicity of Anderson insulator chains with two-body
interactions. A similar problem has been widely studied for systems with
nearest-neighbor density-density interactions, mappable onto the random-field
spin-1/2 XXZ chain. Due to its simplicity, the random-field XXZ chain appears
to be the most natural model for the studies on MBL. However, the interaction
term in the fermionic models mappable onto the XXZ chains is rather specific
in that it has a large projection on (products) of LIOMs of the single-
particle Anderson insulator. Due to this important property, only a tiny part
of the two-body interaction represents a true perturbation to the Anderson
insulators.

One may expect that a generic local two-body interaction also contains other
terms which have no projections on the Anderson LIOMs and thus have greater
impact on the ergodicity of disordered systems. For this reason, in the
present work we have studied two cases of the pair-hopping interactions which
have negligible projections on the Anderson LIOMs.

We found that the results for finite systems at strong disorder $W$ and weak
two-body interaction $V$ qualitatively reproduce the known results for random-
field XXZ chains. However, the regime with exceedingly different model
parameters is most demanding for numerical calculations. For this reason, our
results neither confirm nor contradict the presence of localization for large
$W$ and $V/W\ll 1$. In the opposite case where $V$ and $W$ are of the same
order of magnitude, which was the main focus of this paper, our numerical
results indicate that macroscopic systems are ergodic. The latter observation
follows from our studies of the average level spacing ratios as well as the
fluctuations of the diagonal matrix elements of local observables.

###### Acknowledgements.

We acknowledge the support by the National Science Centre, Poland via project
2020/37/B/ST3/00020 (B.K and M.M.), the support by the Slovenian Research
Agency (ARRS), Research Core Funding Grants P1-0044 (L.V. and J.B.) and
N1-0273 (L.V.). The numerical calculation were partly carried out at the
facilities of the Wroclaw Centre for Networking and Supercomputing.

## References

  * (1) A. Pal and D. A. Huse, Many-body localization phase transition, Phys. Rev. B 82, 174411 (2010).
  * (2) O. S. Barišić and P. Prelovšek, Conductivity in a disordered one-dimensional system of interacting fermions, Phys. Rev. B 82, 161106 (2010).
  * (3) D. J. Luitz, N. Laflorencie, and F. Alet, Many-body localization edge in the random-field Heisenberg chain, Phys. Rev. B 91, 081103 (2015).
  * (4) D. J. Luitz, N. Laflorencie, and F. Alet, Extended slow dynamical regime prefiguring the many-body localization transition, Phys. Rev. B 93, 060201(R) (2016).
  * (5) E. J. Torres-Herrera and L. F. Santos, Dynamics at the many-body localization transition, Phys. Rev. B 92, 014208 (2015).
  * (6) S. Bera, H. Schomerus, F. Heidrich-Meisner, and J. H. Bardarson, Many-body localization characterized from a one-particle perspective, Phys. Rev. Lett. 115, 046603 (2015).
  * (7) J. Hauschild, F. Heidrich-Meisner, and F. Pollmann, Domain-wall melting as a probe of many-body localization, Physical Review B 94, 161109 (2016).
  * (8) T. Devakul and R. R. P. Singh, Early breakdown of area-law entanglement at the many-body delocalization transition, Phys. Rev. Lett. 115, 187201 (2015).
  * (9) C. L. Bertrand and A. M. García-García, Anomalous Thouless energy and critical statistics on the metallic side of the many-body localization transition, Phys. Rev. B 94, 144201 (2016).
  * (10) E. V. H. Doggen, F. Schindler, K. S. Tikhonov, A. D. Mirlin, T. Neupert, D. G. Polyakov, and I. V. Gornyi, Many-body localization and delocalization in large quantum chains, Phys. Rev. B 98, 174202 (2018).
  * (11) J. Šuntajs, J. Bonča, T. Prosen, and L. Vidmar, Quantum chaos challenges many-body localization, Phys. Rev. E 102, 062144 (2020).
  * (12) J. Šuntajs, J. Bonča, T. Prosen, and L. Vidmar, Ergodicity breaking transition in finite disordered spin chains, Phys. Rev. B 102, 064207 (2020).
  * (13) P. Sierant, M. Lewenstein, and J. Zakrzewski, Polynomially filtered exact diagonalization approach to many-body localization, Phys. Rev. Lett. 125, 156601 (2020).
  * (14) D. Basko, I. Aleiner, and B. Altshuler, Metal–insulator transition in a weakly interacting many-electron system with localized single-particle states, Ann. Phys. 321, 1126 (2006).
  * (15) V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature, Phys. Rev. B 75, 155111 (2007).
  * (16) M. Žnidarič, T. Prosen, and P. Prelovšek, Many-body localization in the Heisenberg XXZ magnet in a random field, Phys. Rev. B 77, 064426 (2008).
  * (17) M. Serbyn, Z. Papić, and D. A. Abanin, Universal slow growth of entanglement in interacting strongly disordered systems, Phys. Rev. Lett. 110, 260601 (2013).
  * (18) M. Mierzejewski, J. Herbrych, and P. Prelovšek, Universal dynamics of density correlations at the transition to the many-body localized state, Phys. Rev. B 94, 224207 (2016).
  * (19) H. P. Lüschen, P. Bordia, S. Scherg, F. Alet, E. Altman, U. Schneider, and I. Bloch, Observation of slow dynamics near the many-body localization transition in one-dimensional quasiperiodic systems, Phys. Rev. Lett. 119, 260401 (2017).
  * (20) P. Bordia, H. Lüschen, S. Scherg, S. Gopalakrishnan, M. Knap, U. Schneider, and I. Bloch, Probing slow relaxation and many-body localization in two-dimensional quasiperiodic systems, Phys. Rev. X 7, 041047 (2017).
  * (21) M. Serbyn, Z. Papić, and D. A. Abanin, Thouless energy and multifractality across the many-body localization transition, Phys. Rev. B 96, 104201 (2017).
  * (22) F. Weiner, F. Evers, and S. Bera, Slow dynamics and strong finite-size effects in many-body localization with random and quasiperiodic potentials, Phys. Rev. B 100, 104204 (2019).
  * (23) T. Chanda, P. Sierant, and J. Zakrzewski, Time dynamics with matrix product states: Many-body localization transition of large systems revisited, Phys. Rev. B 101, 035148 (2020).
  * (24) T. C. Berkelbach and D. R. Reichman, Conductivity of disordered quantum lattice models at infinite temperature: Many-body localization, Phys. Rev. B 81, 224429 (2010).
  * (25) Y. Bar Lev, G. Cohen, and D. R. Reichman, Absence of diffusion in an interacting system of spinless fermions on a one-dimensional disordered lattice, Phys. Rev. Lett. 114, 100601 (2015).
  * (26) O. S. Barišić, J. Kokalj, I. Balog, and P. Prelovšek, Dynamical conductivity and its fluctuations along the crossover to many-body localization, Phys. Rev. B 94, 045126 (2016).
  * (27) R. Steinigeweg, J. Herbrych, F. Pollmann, and W. Brenig, Typicality approach to the optical conductivity in thermal and many-body localized phases, Phys. Rev. B 94, 180401 (2016).
  * (28) S. Bera, G. De Tomasi, F. Weiner, and F. Evers, Density propagator for many-body localization: Finite-size effects, transient subdiffusion, and exponential decay, Phys. Rev. Lett. 118, 196801 (2017).
  * (29) P. Prelovšek, M. Mierzejewski, O. S. Barišić, and J. Herbrych, Density correlations and transport in models of many‐body localization, Annalen der Physik 529, 1600362 (2016).
  * (30) J. H. Bardarson, F. Pollmann, and J. E. Moore, Unbounded growth of entanglement in models of many-body localization, Phys. Rev. Lett. 109, 017202 (2012).
  * (31) J. A. Kjäll, J. H. Bardarson, and F. Pollmann, Many-body localization in a disordered quantum Ising chain, Phys. Rev. Lett. 113, 107204 (2014).
  * (32) E. Baygan, S. P. Lim, and D. N. Sheng, Many-body localization and mobility edge in a disordered spin - 1/2 Heisenberg ladder, Phys. Rev. B 92, 195153 (2015).
  * (33) F. Pietracaprina, G. Parisi, A. Mariano, S. Pascazio, and A. Scardicchio, Entanglement critical length at the many-body localization transition, J. Stat. Mech. (2017), 113102 (2017).
  * (34) P. Sierant, D. Delande, and J. Zakrzewski, Many-body localization due to random interactions, Phys. Rev. A 95, 021601 (2017).
  * (35) P. Sierant and J. Zakrzewski, Level statistics across the many-body localization transition, Phys. Rev. B 99, 104205 (2019).
  * (36) M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Evidence for unbounded growth of the number entropy in many-body localized phases, Phys. Rev. Lett. 124, 243601 (2020).
  * (37) D. Sels and A. Polkovnikov, Dynamical obstruction to localization in a disordered spin chain, Phys. Rev. E 104, 054105 (2021).
  * (38) M. Kiefer-Emmanouilidis, R. Unanyan, M. Fleischhauer, and J. Sirker, Slow delocalization of particles in many-body localized phases, Phys. Rev. B 103, 024203 (2021).
  * (39) T. LeBlond, D. Sels, A. Polkovnikov, and M. Rigol, Universality in the onset of quantum chaos in many-body systems, Phys. Rev. B 104, L201117 (2021).
  * (40) L. Vidmar, B. Krajewski, J. Bonča, and M. Mierzejewski, Phenomenology of spectral functions in disordered spin chains at infinite temperature, Phys. Rev. Lett. 127, 230603 (2021).
  * (41) D. Sels and A. Polkovnikov, Thermalization of dilute impurities in one-dimensional spin chains, Phys. Rev. X 13, 011041 (2023).
  * (42) D. Sels, Bath-induced delocalization in interacting disordered spin chains, Phys. Rev. B 106, L020202 (2022).
  * (43) F. Evers and S. Bera, The internal clock of many-body (de-)localization, arXiv:2302.11384 (2023).
  * (44) P. Sierant and J. Zakrzewski, Challenges to observation of many-body localization, Phys. Rev. B 105, 224203 (2022).
  * (45) A. Morningstar, L. Colmenarez, V. Khemani, D. J. Luitz, and D. A. Huse, Avalanches and many-body resonances in many-body localized systems, Phys. Rev. B 105, 174205 (2022).
  * (46) D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Colloquium: Many-body localization, thermalization, and entanglement, Rev. Mod. Phys. 91, 021001 (2019).
  * (47) R. K. Panda, A. Scardicchio, M. Schulz, S. R. Taylor, and M. Žnidarič, Can we study the many-body localisation transition?, EPL (Europhysics Letters) 128, 67003 (2020).
  * (48) D. Abanin, J. Bardarson, G. De Tomasi, S. Gopalakrishnan, V. Khemani, S. Parameswaran, F. Pollmann, A. Potter, M. Serbyn, and R. Vasseur, Distinguishing localization from chaos: Challenges in finite-size systems, Ann. Phys. 427, 168415 (2021).
  * (49) P. J. D. Crowley and A. Chandran, A constructive theory of the numerically accessible many-body localized to thermal crossover, SciPost Phys. 12, 201 (2022).
  * (50) N. Bölter and S. Kehrein, Scrambling and many-body localization in the XXZ chain, Phys. Rev. B 105, 104202 (2022).
  * (51) A. De Luca and A. Scardicchio, Ergodicity breaking in a model showing many-body localization, EPL (Europhysics Letters) 101, 37003 (2013).
  * (52) B. Krajewski, L. Vidmar, J. Bonča, and M. Mierzejewski, Restoring ergodicity in a strongly disordered interacting chain, Phys. Rev. Lett. 129, 260601 (2022).
  * (53) J. M. Deutsch, Quantum statistical mechanics in a closed system, Phys. Rev. A 43, 2046 (1991).
  * (54) M. Srednicki, Chaos and quantum thermalization, Phys. Rev. E 50, 888 (1994).
  * (55) M. Rigol, V. Dunjko, and M. Olshanii, Thermalization and its mechanism for generic isolated quantum systems, Nature (London) 452, 854 (2008).
  * (56) L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics, Adv. Phys. 65, 239 (2016).
  * (57) Y. Y. Atas, E. Bogomolny, O. Giraud, and G. Roux, Distribution of the ratio of consecutive level spacings in random matrix ensembles, Phys. Rev. Lett. 110, 084101 (2013).
  * (58) V. Khemani, S. P. Lim, D. N. Sheng, and D. A. Huse, Critical properties of the many-body localization transition, Phys. Rev. X 7, 021013 (2017).
  * (59) P. Sierant, M. Lewenstein, and J. Zakrzewski, Polynomially filtered exact diagonalization approach to many-body localization, Phys. Rev. Lett. 125, 156601 (2020).
  * (60) F. Andraschko, T. Enss, and J. Sirker, Purification and many-body localization in cold atomic gases, Phys. Rev. Lett. 113, 217201 (2014).
  * (61) H. Kim, T. N. Ikeda, and D. A. Huse, Testing whether all eigenstates obey the eigenstate thermalization hypothesis, Phys. Rev. E 90, 052105 (2014).
  * (62) R. Mondaini, K. R. Fratus, M. Srednicki, and M. Rigol, Eigenstate thermalization in the two-dimensional transverse field Ising model, Phys. Rev. E 93, 032104 (2016).
  * (63) D. Jansen, J. Stolpp, L. Vidmar, and F. Heidrich-Meisner, Eigenstate thermalization and quantum chaos in the Holstein polaron model, Phys. Rev. B 99, 155130 (2019).
  * (64) M. Mierzejewski, M. Środa, J. Herbrych, and P. Prelovšek, Resistivity and its fluctuations in disordered many-body systems: From chains to planes, Phys. Rev. B 102, 161111 (2020).
  * (65) B. Krajewski, M. Mierzejewski, and J. Bonča, Modeling sample-to-sample fluctuations of the gap ratio in finite disordered spin chains, Phys. Rev. B 106, 014201 (2022).