#  Preferred basis of states derived from eigenstate thermalization hypothesis

Hua Yan  Department of Modern Physics, University of Science and Technology of
China, Hefei 230026, China  CAS Key Laboratory of Microscale Magnetic
Resonance, University of Science and Technology of China, Hefei 230026, China
Jiaozi Wang  Department of Modern Physics, University of Science and
Technology of China, Hefei 230026, China  CAS Key Laboratory of Microscale
Magnetic Resonance, University of Science and Technology of China, Hefei
230026, China  Department of Physics, University of Osnabrück, D-49076
Osnabrück, Germany  Wen-ge Wang<EMAIL_ADDRESS>Department of Modern
Physics, University of Science and Technology of China, Hefei 230026, China
CAS Key Laboratory of Microscale Magnetic Resonance, University of Science and
Technology of China, Hefei 230026, China

###### Abstract

We study the long-time average of the reduced density matrix (RDM) of a
central system that is locally coupled to a large environment, under a
Schrödinger evolution of the total system. We consider a class of interaction
Hamiltonian, whose environmental part satisfies the so-called eigenstate
thermalization hypothesis ansatz with a constant diagonal part in the energy
region concerned. Relations among elements of the averaged RDM are derived.
When steady states of the central system exist, these relations imply the
existence of a preferred basis, which is given by the eigenbasis of a
renormalized self-Hamiltonian that includes certain averaged impact of the
system-environment interaction. Numerical simulations performed for a qubit
coupled to a defect Ising chain confirm the analytical predictions.

##  I Introduction

Properties of small open quantum systems, which are coupled to large quantum
environments, have attracted significant attention and been studied
extensively in recent decades in various fields of physics [1, 2, 3, 4, 5].
Such a central system is described by its reduced density matrix (RDM) and may
approach a steady state in many situations. For example, it is now well known
that the phenomenon of decoherence, due to interactions with huge quantum
environments, may happen in such a way that a RDM becomes approximately
diagonal on a so-called preferred (pointer) basis of states (PBS) [6, 7, 8, 9,
10, 11]. Under pure-dephasing interactions, decoherence has been studied well,
with PBS given by eigenbases of self-Hamiltonians [12, 13, 14]. However, under
strong interactions and complex environments, with the self-Hamiltonians
negligible, PBS may be given by eigenbases of the interaction Hamiltonians [6,
12].

The situation is much more complicated with a generic dissipative interaction,
whose Hamiltonian is not commutable with the central system’s self-
Hamiltonian, due to the interplay of decoherence and relaxation. In this
generic case, knowledge about PBS is still far from being complete. Under a
sufficiently weak interaction and by a first-order perturbation theory, it was
found that the system’s eigenbasis is approximately a PBS under a quantum
chaotic environment [15]. When the total system’s eigenfunctions possess
certain special randomness, a PBS (if existing) is given by the eigenbasis of
a renormalized self-Hamiltonian [16]. These results are in agreement with a
generic expectation for Markovian processes described by Lindblad master
equations, as exemplified in solvable models [2, 3, 13]. While, when non-
Markovian effects due to dynamics of the total system are taken into account,
nonnegligible off-diagonal elements of RDM have been found at long times in
various spin-boson models on the spin’s energy basis [17, 18, 19, 20, 21].

In this paper, we go further in the study of properties of steady states of
small open systems, by directly computing their long-time averaged RDM under
overall Schrödinger evolutions of total systems. A key point of our approach
is to consider those environments, for which the environmental parts of the
interaction Hamiltonians satisfy the so-called eigenstate thermalization
hypothesis (ETH) ansatz [22, 23, 24, 25, 26, 27] and their diagonal elements
in the ansatz may be treated as constants within the energy regions of
relevance. We are to derive $\frac{1}{2}(m-1)(m+2)$ relations among elements
of such an averaged RDM for a central system with a number of $m$ levels. When
steady states exist, these relations imply that the central system should have
a PBS, which is given by the eigenbasis of a renormalized self-Hamiltonian
that includes certain impact of the system-environment interaction.

The paper is organized as follows. In Sec. II, we specify the systems to be
studied. In Sec. III, we derive the above-mentioned relations. Some further
discussions are given in Sec.IV. Numerical simulations are presented in Sec.
V, to illustrate validity of the analytical predictions for a qubit as the
central system and a defect Ising chain as the environment. Finally,
conclusions and discussions are given in Sec. VI.

##  II Setup

In this section, we discuss basic properties of the Hamiltonians of the
systems to be studied. We use $S$ to denote the central system and use
$\mathcal{E}$ to denote its (large) environment which consists of $N$
particles ($N\gg 1$). Hilbert spaces of $S$ and $\mathcal{E}$ are denoted by
$\mathcal{H}^{S}$ and $\mathcal{H}^{\mathcal{E}}$, respectively, with
dimensions $m$ and $d_{\mathcal{E}}$. The value of $m$ is required to be much
smaller than the number of environmental levels that are of relevance
effectively to the time evolution.

The Hamiltonian of the total system is written as

$H=H^{S}+H^{I}+H^{\mathcal{E}},$ (1)

where $H^{S}$ and $H^{\mathcal{E}}$ are self-Hamiltonians of $S$ and
$\mathcal{E}$, respectively, which are obtained in the weak coupling limit,
and $H^{I}$ represents a local interaction Hamiltonian. Normalized eigenstates
of the total system are denoted by $|n\rangle$ with energies $E_{n}$ in the
increasing-energy order,

$\displaystyle H|n\rangle=E_{n}|n\rangle.$ (2)

Normalized eigenstates of $H^{S}$ and of $H^{\mathcal{E}}$ are denoted by
$|\alpha\rangle$ and $|i\rangle$, respectively, with labels $\alpha$ and $i$
as positive integers starting from $1$. The corresponding eigenenergies are
denoted by $e^{S}_{\alpha}$ and $e_{i}$, respectively, both in the increasing-
energy order,

$\displaystyle H^{S}|\alpha\rangle=e^{S}_{\alpha}|\alpha\rangle,$ (3a)
$\displaystyle H^{\mathcal{E}}|i\rangle=e_{i}|i\rangle,$ (3b)

where for brevity we have omitted a superscript $\mathcal{E}$ for $e_{i}$. We
use $\Delta_{S}$ to indicate the energy scope the central system $S$:

$\displaystyle\Delta_{S}:=e^{S}_{m}-e^{S}_{1}.$ (4)

We use $H^{0}$ to indicate the uncoupled Hamiltonian,

$\displaystyle H^{0}=H^{S}+H^{\mathcal{E}}.$ (5)

Its eigenstates are written as $|\alpha\rangle|i\rangle$, in short, $|\alpha
i\rangle$, satisfying $H^{0}|\alpha i\rangle=E_{\alpha i}|\alpha i\rangle$,
where $E_{\alpha i}=e^{S}_{\alpha}+e_{i}$. The expansion of a state
$|n\rangle$ on the basis given by $|\alpha i\rangle$ is written as

$|n\rangle=\sum_{\alpha i}C_{\alpha i}^{n}|\alpha i\rangle,$ (6)

with expansion coefficients $C_{\alpha i}^{n}$. For the simplicity in
discussion, we consider a product form of $H^{I}$, 111Generalization to a
generic local interaction Hamiltonian will be briefly discussed in Sec. IV.2.

$H^{I}=\lambda H^{IS}\otimes H^{I\mathcal{E}},$ (7)

where $H^{IS}$ and $H^{I\mathcal{E}}$ are Hermitian operators acting on the
two spaces of $\mathcal{H}^{S}$ and $\mathcal{H}^{\mathcal{E}}$, respectively,
and $\lambda$ is a parameter for characterizing the coupling strength.
Elements of $H^{IS}$ and $H^{I\mathcal{E}}$ on $|\alpha\rangle$ and
$|i\rangle$ are written as

$\displaystyle H^{IS}_{\alpha\beta}=\langle\alpha|H^{IS}|\beta\rangle,$ (8a)
$\displaystyle H^{I\mathcal{E}}_{ij}=\langle i|H^{I\mathcal{E}}|j\rangle.$
(8b)

To describe locality of the interaction, we further divide the environment
$\mathcal{E}$ into a small part denoted by $\mathcal{E}_{1}$ and a large part
denoted by $\mathcal{E}_{2}$, such that the system $S$ is coupled to
$\mathcal{E}_{1}$ only. Then, $H^{I\mathcal{E}}$ is written as

$\displaystyle H^{I\mathcal{E}}=H^{I\mathcal{E}_{1}}\otimes
I^{\mathcal{E}_{2}},$ (9)

where $H^{I\mathcal{E}_{1}}$ is an operator that acts on the Hilbert space of
$\mathcal{E}_{1}$ and $I^{\mathcal{E}_{2}}$ indicates the identity operator on
the Hilbert space of $\mathcal{E}_{2}$. 222As a local operator,
$H^{I\mathcal{E}_{1}}$ does not change with the environmental particle number
$N$.

Although the exact condition under which the ETH ansatz proposed in Ref.[24]
is applicable is still unclear, it is usually expected valid at least for
local operators of many-body quantum chaotic systems [26, 27]. Here, we assume
that this ansatz is applicable to the operator $H^{I\mathcal{E}}$. According
to this hypothesis, (1) diagonal elements $H^{I\mathcal{E}}_{ii}$ on average
vary slowly with the eigenenergy $e_{i}$; (2) fluctuations of
$H^{I\mathcal{E}}_{ii}$ possess certain random feature and are very small,
scaling as $e^{-S(e)/2}$, where $S(e)$ is proportional to the particle number
$N$ of $\mathcal{E}$ and is related to the micro-canonical entropy in a
semiclassical treatment; and (3) off-diagonal elements $H^{I\mathcal{E}}_{ij}$
with $i\neq j$ behave in a way similar to fluctuations of
$H^{I\mathcal{E}}_{ii}$ [22, 23, 25, 24, 26, 28]. These predictions are
written in the following concise form, usually referred to as the ETH ansatz,

$H^{I\mathcal{E}}_{ij}={h}(e)\delta_{ij}+e^{-S(e)/2}g(e,\omega)R_{ij},$ (10)

where $e=(e_{i}+e_{j})/2$, $\omega=e_{j}-e_{i}$, ${h}(e)$ is a slowly varying
function of $e$, $g(e,\omega)$ is some smooth function, and $R_{ij}$ indicate
random variables with a normal distribution (zero mean and unit variance).
333Certain correlations among $R_{ij}$ have been observed numerically in some
chaotic systems [29], but, we do not discuss this possibility in this paper.

Analytical expressions of the functions $h(e)$ and $g(e,\omega)$ are still
lacking. Numerically, three regimes have been observed for $|g(e,\omega)|$
with respect to the order of per-site energy denoted by $\xi$, provided that
$e$ lies in the central region of the spectrum [26, 5]. That is, for
$\omega\ll\xi$, it shows a plateau with a height proportional to $N^{1/2}$ and
a width proportional to $N^{-2}$ [25, 30]; for large $\omega\gg\xi$, it decays
exponentially; and, for $\omega\sim\xi$, it is proportional to $\omega^{-1/2}$
in diffusive one-dimensional systems [31, 32, 33, 34].

For the simplicity in discussion, we set the initial state of the total system
at a time $t=0$ as a pure state with a product form, 444Discussions to be
given below may be generalized, in a straightforward way, to a generic initial
state written as
$|\Psi(0)\rangle=\sum_{\alpha}c_{0\alpha}|\alpha\rangle\otimes|\mathcal{E}_{0\alpha}\rangle$,
if all the environmental states $|\mathcal{E}_{0\alpha}\rangle$ lie in the
same energy shell $\Gamma^{\mathcal{E}}_{0}$ in Eq.(14).  that is,

$|\Psi(0)\rangle=|\phi_{S}\rangle\otimes|\mathcal{E}_{0}\rangle.$ (11)

Here, $|\phi_{S}\rangle$ indicates an arbitrary normalized state of the
central system $S$, written as

$\displaystyle|\phi_{S}\rangle=\sum_{\alpha}c_{0\alpha}|\alpha\rangle;$ (12)

and $|\mathcal{E}_{0}\rangle$ is an arbitrary environmental state that lies
within an energy shell denoted by $\Gamma^{\cal E}_{0}$,

$|\mathcal{E}_{0}\rangle=\sum_{e_{i}\in\Gamma^{\cal E}_{0}}c_{0i}|i\rangle.$
(13)

The energy shell $\Gamma^{\cal E}_{0}$ is centered at an energy $e_{0}$ and
has a width $\delta e_{0}$, 555 See Sec.III.4 for a discussion about
restriction to the width $\delta e_{0}$.  namely,

$\displaystyle\Gamma^{\cal E}_{0}=[e_{0}-\delta e_{0}/2,e_{0}+\delta
e_{0}/2].$ (14)

##  III Main result

In this section, we derive the main result of this paper, as relations among
elements of the long-time averaged RDM. Specifically, we give some formal
discussions in Sec.III.1, then, in Sec.III.2, derive an upper bound to the
environmental energy region, which is of relevance effectively to the wave
function at all times. The main result is derived in Sec.III.3 and properties
of a main condition used in it are discussed in Sec.III.4.

###  III.1 Preliminary discussions

The total system undergoes a Schrödinger evolution,

$\displaystyle|\Psi(t)\rangle=e^{-iHt/\hbar}|\Psi(0)\rangle.$ (15)

We write $|\Psi(t)\rangle$ in the following expansion with respect to the
central system’s states $|\alpha\rangle$,

$\displaystyle|\Psi(t)\rangle=\sum_{\alpha=1}^{m}|\alpha\rangle|\mathcal{E}_{\alpha}(t)\rangle,$
(16)

and call $|\mathcal{E}_{\alpha}(t)\rangle$ the _environmental branches_ of
$|\Psi(t)\rangle$. These branches, as vectors in the environmental Hilbert
space, are written as

$\displaystyle|\mathcal{E}_{\alpha}(t)\rangle=\langle\alpha|\Psi(t)\rangle,$
(17)

and are usually not normalized. Under the initial condition in Eq.(11), it is
direct to find that

$|\mathcal{E}_{\alpha}(t)\rangle=\sum_{\beta}c_{0\beta}\langle\alpha|e^{-iHt/\hbar}|\beta\rangle|\mathcal{E}_{0}\rangle,$
(18)

and

$i\hbar\frac{d}{dt}|\mathcal{E}_{\alpha}(t)\rangle=H_{\alpha\alpha}|\mathcal{E}_{\alpha}(t)\rangle+\sum_{\beta\neq\alpha}H_{\alpha\beta}|\mathcal{E}_{\beta}(t)\rangle,$
(19)

where $H_{\alpha\beta}$ indicate operators that act on the Hilbert space of
the environment, as defined below,

$\displaystyle H_{\alpha\beta}:=\langle\alpha|H|\beta\rangle.$ (20)

By definition, the RDM of the system $S$, denoted by $\rho^{S}(t)$, is given
by $\rho^{S}(t)=\text{Tr}_{\mathcal{E}}\rho(t)$, where
$\rho(t)=|\Psi(t)\rangle\langle\Psi(t)|$. It is easy to check that elements of
the RDM on the basis $\\{|\alpha\rangle\\}$, namely
$\rho_{\alpha\beta}^{S}(t)=\langle\alpha|\rho^{S}(t)|\beta\rangle$, have the
following expression,

$\rho_{\alpha\beta}^{S}(t)=\langle\mathcal{E}_{\beta}(t)|\mathcal{E}_{\alpha}(t)\rangle.$
(21)

Making use of Eq.(19), after some deviation, one finds that the elements
$\rho_{\alpha\beta}^{S}(t)$ satisfy the following equation (see Appendix A),

$\displaystyle
i\hbar\frac{d\rho_{\alpha\beta}^{S}(t)}{dt}=W_{\alpha\beta}^{(1)}+\lambda
W_{\alpha\beta}^{(2)},$ (22)

where

$\displaystyle W_{\alpha\beta}^{(1)}$ $\displaystyle=$
$\displaystyle(e_{\alpha}^{S}-e_{\beta}^{S})\rho_{\alpha\beta}^{S}(t),$ (23a)
$\displaystyle W_{\alpha\beta}^{(2)}$ $\displaystyle=$
$\displaystyle\sum_{\gamma=1}^{m}H_{\alpha\gamma}^{IS}F_{\beta\gamma}(t)-\sum_{\gamma=1}^{m}H^{IS}_{\gamma\beta}F_{\gamma\alpha}(t).$
(23b)

Here, $F_{\alpha\beta}(t)$ indicate $c$-number quantities defined below,

$F_{\alpha\beta}(t):=\langle\mathcal{E}_{\alpha}(t)|H^{I\mathcal{E}}|\mathcal{E}_{\beta}(t)\rangle,$
(24)

and, from them, we define the following operator,

$\displaystyle
F(t):=\sum_{\alpha\beta}F_{\alpha\beta}(t)|\alpha\rangle\langle\beta|.$ (25)

It is easy to check that $W_{\alpha\beta}^{(1)}$ and $W_{\alpha\beta}^{(2)}$
have the following concise expressions,

$\displaystyle
W_{\alpha\beta}^{(1)}=\langle\alpha|[H^{S},\rho^{S}(t)]|\beta\rangle,$ (26a)
$\displaystyle
W_{\alpha\beta}^{(2)}=\langle\alpha|[H^{IS},F^{T}(t)]|\beta\rangle,$ (26b)

where ${F}^{T}$ indicates the transposition operator of ${F}$, which is
defined on the eigenbasis of $H^{S}$.

We use an overline to indicate the long-time average of a term. For example,
the long-time average of the RDM is written as $\overline{\rho}^{S}$,

$\overline{\rho}^{S}=\lim_{t\to\infty}\frac{1}{t}\int_{0}^{t}\rho^{S}(t^{\prime})dt^{\prime}.$
(27)

Clearly, in the case that a steady state of the RDM exists, it is given by
$\overline{\rho}^{S}$. Since the elements $\rho_{\alpha\beta}^{S}(t)$ have
bounded values, the long-time average of $d\rho_{\alpha\beta}^{S}(t)/{dt}$
must be zero, i.e., $\overline{{d\rho_{\alpha\beta}^{S}(t)}/{dt}}=0$. Then,
Eq.(22) gives that

$\displaystyle\overline{W}_{\alpha\beta}^{(1)}+\lambda\overline{W}_{\alpha\beta}^{(2)}=0.$
(28)

Substituting the explicit expressions of ${W}_{\alpha\beta}^{(1)}$ and
${W}_{\alpha\beta}^{(2)}$ in Eq.(23) into Eq.(28), one finds the following
formal relation for the long-time averaged RDM:

$\displaystyle(e_{\alpha}^{S}-e_{\beta}^{S})\overline{\rho}_{\alpha\beta}^{S}+\lambda\sum_{\gamma=1}^{m}\Big{[}H_{\alpha\gamma}^{IS}\overline{F}_{\beta\gamma}-H^{IS}_{\gamma\beta}\overline{F}_{\gamma\alpha}\Big{]}=0.$
(29)

It is straightforward to check that a concise form of Eq.(29) is written as

$\displaystyle[H^{S},\overline{\rho}^{S}]+\lambda[H^{IS},\overline{F}^{T}]=0.$
(30)

###  III.2 Effective environmental energy region

In this section, we discuss an environmental energy region, within which all
the branches $|\mathcal{E}_{\alpha}(t)\rangle$ effectively lie for all the
times $t$, and indicate it by $\Gamma^{\mathcal{E}}$. We do not need to find
the smallest one of this type of region. Instead, we consider a region that
has the following simple form,

$\displaystyle\Gamma^{\mathcal{E}}=[e_{0}-\delta e/2,e_{0}+\delta e/2],$ (31)

centered at the initial center $e_{0}$ and with a width $\delta e$.

Below, we derive an expression for $\delta e$, as an upper bound to the width
of the energy region that effectively contains all
$|\mathcal{E}_{\alpha}(t)\rangle$. For this purpose, we need to analyze the
components $\langle i|\mathcal{E}_{\alpha}(t)\rangle$,

$\displaystyle\langle
i|\mathcal{E}_{\alpha}(t)\rangle=\sum_{\beta,n}\sum_{e_{j}\in\Gamma_{0}^{\mathcal{E}}}c_{0\beta}c_{0j}C^{n*}_{\beta
j}C^{n}_{\alpha i}e^{-iE_{n}t},$ (32)

which is directly obtained by making use of Eqs.(11)-(13) and (15)-(17).
Initially, with $e^{-iE_{n}t}=1$ at $t=0$, due to correlations among the terms
of $(C^{n*}_{\beta j}C^{n}_{\alpha i})$ of different indices $n$, which
originate from the completeness of the states $|n\rangle$ as a basis in the
total Hilbert space, nonzero values of the rhs of Eq.(32) are restricted
within the initial energy region $\Gamma^{\mathcal{E}}_{0}$. With increase of
the time $t$, the phases $e^{-iE_{n}t}$ gradually destroy the above-mentioned
correlations and, as a result, the energy region that is effectively occupied
by $|\mathcal{E}_{\alpha}(t)\rangle$ expands. Cutting all the correlations by
taking an absolute value for each summed term on the rhs of Eq.(32), one gets
an upper bound to $|\langle i|\mathcal{E}_{\alpha}(t)\rangle|$:

$\displaystyle|\langle
i|\mathcal{E}_{\alpha}(t)\rangle|\leq\sum_{e_{j}\in\Gamma_{0}^{\mathcal{E}}}\sum_{\beta,n}\left|c_{0\beta}c_{0j}C^{n*}_{\beta
j}C^{n}_{\alpha i}\right|.$ (33)

When using the rhs of Eq.(33) to get an upper bound to the environmental
energy region that effectively contains $|\mathcal{E}_{\alpha}(t)\rangle$,
exact values of the nonzero coefficients $c_{0\beta}$ and $c_{0j}$ are not
important. Hence, we may focus on the values of $C^{n*}_{\beta j}C^{n}_{\alpha
i}$. In particular, for the eigenfunction (EF) $C^{n}_{\alpha i}$ of each
state $|n\rangle$, what is of relevance is its main-body region, within which
the main population lie up to a small error indicated by $\epsilon$.
Energetically, such a main-body region consists of those uncoupled states
$|\alpha i\rangle$, whose energies $E_{\alpha i}$ are around the exact energy
$E_{n}$ within a scope which we indicate by $w^{\epsilon}_{n}$. More exactly,
the set of the indices of these uncoupled states, indicated by
$\Omega^{\epsilon}_{n}$, is written as

$\displaystyle\Omega^{\epsilon}_{n}=\left\\{(\alpha,i):|E_{\alpha
i}-E_{n}|\leq\frac{1}{2}w^{\epsilon}_{n}\right\\}.$ (34)

Then, the main-body region of the EF $C^{n}_{\alpha i}$ satisfies the
following requirement,

$\displaystyle\sum_{(\alpha,i)\in\Omega^{\epsilon}_{n}}\left|C^{n}_{\alpha
i}\right|^{2}\doteq 1-\epsilon\quad(\epsilon\ll 1),$ (35)

where “$\doteq$” means that the set $\Omega^{\epsilon}_{n}$ is chosen such
that the left hand side of Eq.(35) is the closest to its right hand side.

Figure 1: A schematic illustration for the energy region
$\Gamma^{\mathcal{E}}$. Upper panel: the environmental branch that moves to
the farthest left from the initial shell $\Gamma^{\mathcal{E}}_{0}$ due to the
system-environment interaction. Lower panel: the branch that moves to the
farthest right.

For a product $C^{n*}_{\beta j}C^{n}_{\alpha i}$ to give a nonnegligible
contribution to the rhs of Eq.(33), both of the two basis states $|\beta
j\rangle$ and $|\alpha i\rangle$ should lie in the main-body region of
$|n\rangle$. As as result, the expansion from $\Gamma^{\mathcal{E}}_{0}$ to
$\Gamma^{\mathcal{E}}$ should be influenced mainly by two factors: widths of
main-body regions of the EFs and the central system’s energy differences. We
use $w^{\epsilon}_{\rm max}$ to indicate the maximum value of
$w^{\epsilon}_{n}$ for those states $|n\rangle$ that are of relevance to the
time evolution of the initial state. Then, noting that $\Delta_{S}$ in Eq.(4)
gives the maximum of $|e^{S}_{\alpha}-e^{S}_{\beta}|$, we find the following
expression of $\delta e$,

$\displaystyle\delta e=\delta e_{0}+2\Delta_{S}+w^{\epsilon}_{\rm max},$ (36)

as illustrated in Fig.1.

Two remarks are in order: (1) At $\lambda=0$, the real width is just $\delta
e_{0}$, smaller than $\delta e$. (2) When the states $|\alpha i\rangle$ are
sufficiently coupled by the interaction, it is possible for $\delta e$ in
Eq.(36) to be close to the width of the energy region that is really occupied.

###  III.3 Relations among elements of averaged RDM

In this section, we derive the main result, as relations that the elements of
$\overline{\rho}^{S}$ satisfy. To this end, let us expand the environmental
branches $|\mathcal{E}_{\alpha}(t)\rangle$ as,

$\displaystyle|\mathcal{E}_{\alpha}(t)\rangle=\sum_{i}f_{\alpha
i}(t)|i\rangle$ (37)

with expansion coefficients $f_{\alpha i}(t)$. Substituting Eq.(37) into
Eqs.(21) and (24), taking the long-time average, and making use of the fact
that all the environmental branches effectively lie within the energy region
$\Gamma^{\mathcal{E}}$, one finds the following expressions of
$\overline{\rho}_{\alpha\beta}^{S}$ and $\overline{F}_{\beta\alpha}$,

$\displaystyle\overline{\rho}_{\alpha\beta}^{S}=\sum_{i}\overline{f_{\beta
i}^{*}f_{\alpha i}}\simeq\sum_{e_{i}\in\Gamma^{\mathcal{E}}}\overline{f_{\beta
i}^{*}f_{\alpha i}},$ (38)
$\displaystyle\overline{F}_{\beta\alpha}=\sum_{i,j}\overline{f_{\beta
j}^{*}f_{\alpha
i}}H^{I\mathcal{E}}_{ji}\simeq\sum_{e_{i},e_{j}\in\Gamma^{\mathcal{E}}}\overline{f_{\beta
j}^{*}f_{\alpha i}}H^{I\mathcal{E}}_{ji}.$ (39)

Substituting the ETH ansatz (10) into Eq.(39), one gets that

$\displaystyle\overline{F}_{\beta\alpha}\simeq\sum_{e_{i}\in{\Gamma}^{\mathcal{E}}}h(e_{i})\overline{f_{\beta
i}^{*}f_{\alpha i}}+\Delta_{\alpha\beta},$ (40)

where $\Delta_{\alpha\beta}$ indicates a fluctuation term, given by

$\displaystyle\Delta_{\alpha\beta}=\sum_{i,j}\overline{f_{\beta
i}^{*}f_{\alpha j}}e^{-S(e)/2}g(e,\omega)R_{ij}.$ (41)

Generically, the two operators $\overline{F}$ and $\overline{\rho}^{S}$ do not
have a simple relationship. One key observation made here is that they may
possess a simple relationship, if the function $h(e)$ is approximately a
nonzero constant within the energy shell $\Gamma^{\mathcal{E}}$. More exactly,
the condition is that

$\displaystyle\left|\frac{1}{h_{0}}\Delta
h\right|\leq\epsilon_{h}\quad\text{with $\epsilon_{h}\ll 1$},$ (42)

where $\epsilon_{h}$ is a parameter much smaller than $1$, $h_{0}\equiv
h(e_{0})$ with $h_{0}\neq 0$ (even in the limit of $N\to\infty$), and $\Delta
h$ indicates the maximum difference between $h_{0}$ and $h(e_{i})$ within
$\Gamma^{\mathcal{E}}$:

$\displaystyle\Delta h=\max_{e_{i}\in{\Gamma}^{\mathcal{E}}}|h(e_{i})-h_{0}|.$
(43)

To show the above-mentioned relationship, we note that, when the condition in
Eq.(42) is satisfied, Eqs.(38) and (40) imply that

$\displaystyle\overline{F}_{\beta\alpha}\simeq
h_{0}\overline{\rho}^{S}_{\alpha\beta}+\Delta_{\alpha\beta},$ (44)

or, equivalently,

$\overline{F}^{T}\simeq h_{0}\overline{\rho}^{S}+\Delta,$ (45)

where $\Delta$ is a fluctuation operator,

$\displaystyle\Delta:=\sum_{\alpha,\beta}\Delta_{\alpha\beta}|\alpha\rangle\langle\beta|.$
(46)

From the rhs of Eq.(41), one sees three factors that influence the
$N$-dependence of $\Delta$. The first one is the exponential decay of
$e^{-S(e)/2}$, with $S(e)\sim N$. The second factor is given by the unknown
ETH-ansatz function $g(e,\omega)$, for which numerical simulations show a
polynomial increase of $N^{\gamma}$, with $\gamma=1/2$ in some (diffusive one-
dimensional) systems [35, 26]. The third factor lies in the summation over the
indices $i$ and $j$ and the long-time average term $\overline{f_{\beta
i}^{*}f_{\alpha j}}$. As shown in Appendix B, due to the randomness of
$R_{ij}$, contribution from the third factor is negligible compared with the
first factor. (See Eq.(97) for an upper bound to the norm of $\Delta$.)
Therefore, the $N$-scaling behavior of the fluctuation operator $\Delta$ is
dominated by the exponential decay term $e^{-S(e)/2}$. Due to this exponential
decay, as well as the fact that $h_{0}\neq 0$ in the limit of large $N$ and
that each RDM has a unit trace, one gets that

$\displaystyle\|\Delta\|\ll\|h_{0}\overline{\rho}^{S}\|\quad\text{at
sufficiently large $N$}.$ (47)

Then, Eq.(45) gives the following relation between the two operators of
$\overline{F}$ and $\overline{\rho}$:

$\displaystyle\overline{F}^{T}\simeq h_{0}\overline{\rho}^{S}\quad\text{at
sufficiently large $N$.}$ (48)

The main result of this paper is obtained by substituting Eq.(48) into
Eq.(30),

$\displaystyle[\widetilde{H}^{S},\overline{\rho}^{S}]\simeq 0,$ (49)

which holds under the condition of Eq.(42) and at sufficiently large $N$.
Here, $\widetilde{H}^{S}$ is a renormalized self-Hamiltonian of the central
system, defined by

$\displaystyle\widetilde{H}^{S}:=H^{S}+\lambda h_{0}H^{IS},$ (50)

which includes certain averaged impact of the system-environment interaction.
From Eq.(49), one sees that, if a PBS exists, it should be given by the
eigenbasis of the renormalized self-Hamiltonian $\widetilde{H}^{S}$. Writing
Eq.(49) explicitly, one gets that

$\displaystyle(e_{\alpha}^{S}-e_{\beta}^{S})\overline{\rho}^{S}_{\alpha\beta}+\lambda
h_{0}\sum_{\gamma}[H^{IS}_{\alpha\gamma}\overline{\rho}^{S}_{\gamma\beta}-\overline{\rho}^{S}_{\alpha\gamma}H^{IS}_{\gamma\beta}]\simeq
0.$ (51)

This gives $m(m-1)/2$ relations among elements of the averaged RDM for
$\alpha\neq\beta$ and $(m-1)$ relations for $\alpha=\beta$.

As an illustration of the above result, let us consider a nondegenerate two-
level system (TLS), with $e_{2}^{S}\neq e_{1}^{S}$. From Eq.(51) with
$\alpha\neq\beta$, one gets that

$\displaystyle\overline{\rho}^{S}_{12}\simeq\frac{\lambda\eta_{r}h_{0}}{1-\lambda\eta_{d}h_{0}}(\overline{\rho}^{S}_{22}-\overline{\rho}^{S}_{11}),$
(52)

where

$\displaystyle\eta_{d}=\frac{H_{11}^{IS}-H_{22}^{IS}}{e_{2}^{S}-e_{1}^{S}},\quad\eta_{r}=\frac{H_{12}^{IS}}{e_{2}^{S}-e_{1}^{S}}.$
(53)

The quantity $\eta_{d}$ gives a relative measure for the strength of
dephasing, while, $\eta_{r}$ gives a relative measure for the strength of
relaxation (dissipation). Meanwhile, in the case of $\alpha=\beta$, one gets
that

$H^{IS}_{12}\overline{\rho}^{S}_{21}-\overline{\rho}^{S}_{12}H^{IS}_{21}\simeq
0,$ (54)

which implies approximate realness of the product
$H^{IS}_{12}\overline{\rho}^{S}_{21}$.

###  III.4 $N$\- and $\lambda$-relevance to the condition (42)

In this section, we discuss relevance of the particle number $N$ to Eq.(42), a
main prerequisite for the above-derived main result, as well as relevance of
the interaction strength $\lambda$. 666 One may note that Eq.(42) is always
satisfied, in the case that EFs of the quantum chaotic environment may be
effectively described by the random matrix theory (RMT). In fact, in this
case, $h(e)$ is a constant, given by
$h(e)=\text{tr}(H^{I\mathcal{E}})/d_{\mathcal{E}}$ [26].  Basically, Eq.(42)
requires that the environmental energy shell $\Gamma^{\mathcal{E}}$ should be
“sufficiently narrow”, such that the function $h(e)$ may be approximately
taken as a constant within it, compared with its nonzero central value
$h_{0}$. Below, we give a detailed discussion of the exact meaning of “being
sufficiently narrow”.

####  III.4.1 Relevance of the particle number N

Relevance of $N$ to Eq.(42) comes mainly from two aspects: the width $\delta
e$ of $\Gamma^{\mathcal{E}}$ in Eq.(36) and the ETH-ansatz function $h(e)$.
The width $\delta e=\delta e_{0}+2\Delta_{S}+w^{\epsilon}_{\rm max}$ contains
three terms. Clearly, $\Delta_{S}$, the central system’s energy scope, is
$N$-independent. The $N$-dependence of $\delta e_{0}$ is usually determined
according to the problem at hand, particularly, to quantities of final
interest; e.g., it may be taken as a constant, or as some polynomial function
of $N$.

The situation with $w^{\epsilon}_{\rm max}$, the maximum width of relevant EFs
of the total system on the uncoupled energy basis, is more complicated. In
fact, presently, still not much is known analytically about widths of the EFs.
It seems reasonable to assume that $w^{\epsilon}_{\rm max}\sim N^{\mu}$ with
some parameter $\mu$ the value of which may be model-dependent. By a first-
order perturbation-theory treatment to long tails of EFs in certain model, it
was found that $\mu<0$ [36]; while, a study of higher-order contributions is
still under investigation [37] by making use of a semiperturbative theory [38,
39, 40, 41].

The ETH ansatz does not assume any specific form of the function $h(e)$.
According to numerical simulations with the help of some analytical analysis
[42, 33, 26, 43], $h(e)$ was found approximately a function of per-site
energy,

$\displaystyle h(e_{i})\approx\widetilde{h}(e_{i}/N),$ (55)

where $\widetilde{h}(x)$ is some smooth function of $x$, independent of $N$.
Then, Taylor’s expansion gives that

$\displaystyle
h(e_{i})-h(e_{j})=\widetilde{h}^{\prime}(e_{i}/N)\frac{e_{i}-e_{j}}{N}+O_{2}(\frac{e_{i}-e_{j}}{N}),$
(56)

where $\widetilde{h}^{\prime}(x)$ indicates the derivative of
$\widetilde{h}(x)$ and $O_{2}$ represents the second and higher order terms of
the expansion.

To be specific, let us discuss a case, in which the initial width $\delta
e_{0}$ increases slower than $N$ such that

$\displaystyle\lim_{N\to\infty}\frac{\delta e_{0}}{N}=0.$ (57)

This case may be met quite often practically. Note that Eq.(57) does not
really require narrowness of the initial shell $\Gamma^{\mathcal{E}}_{0}$;
e.g., it holds for $\delta e_{0}\sim N^{b}$ with a parameter $b<1$. Then, as
long as $\mu<1$ for $w^{\epsilon}_{\rm max}\sim N^{\mu}$, Eq.(57) implies that

$\displaystyle\lim_{N\to\infty}\frac{\delta e}{N}=0.$ (58)

This implies that the ratio $|e_{i}-e_{j}|/{N}$ should approach zero in the
limit of large $N$ for $e_{i},e_{j}\in\Gamma^{\mathcal{E}}$. If
$\widetilde{h}^{\prime}\neq 0$, then, according to Eq.(56), the difference
$[h(e_{i})-h(e_{j})]$ is approximately given by the first-order term at
sufficiently large $N$. As a consequence, $h(e)$ is approximately a linear
function within $\Gamma^{\mathcal{E}}$ and $\Delta h$ in Eq.(43) is written as

$\displaystyle\Delta h\approx\frac{\delta
e}{2N}\left|\widetilde{h}^{\prime}\left(\frac{e_{0}}{N}\right)\right|.$ (59)

One sees that, as long as
$\left|\widetilde{h}^{\prime}\left({e_{0}}/N\right)\right|$ has a finite upper
bound, $\Delta h/h_{0}\to 0$ in the limit of $N\to\infty$. Otherwise, i.e., if
$\widetilde{h}^{\prime}=0$, one may consider the second-order term (if
nonzero) in Eq.(56) and, following arguments similar to those given above,
reach the same conclusion. Similar arguments also apply, when higher-order
terms dominate. Therefore, for systems with $\mu<1$, under an initial
condition satisfying Eq.(57), the condition (42) is usually fulfilled at
sufficiently large $N$.

####  III.4.2 Relevance of the interaction strength

Among the three terms of $\Delta_{S}$, $\delta e_{0}$, and $w^{\epsilon}_{\rm
max}$ in $\delta e$, only the EF width $w^{\epsilon}_{\rm max}$ depends on the
interaction strength $\lambda$. As is well known, usually, $w^{\epsilon}_{\rm
max}$ increases with increasing $\lambda$, when other parameters in the total
Hamiltonian are fixed. It is reasonable to expect that dependence of
$w^{\epsilon}_{\rm max}$ on the pair of $(N,\lambda)$ may behave in a quite
complicated way. A full understanding of this behavior is beyond the scope of
this investigation. Below, for the sake of clearness in discussion, we usually
consider a fixed value of $N$ when discussing influence of $\lambda$.

To study influence of the interaction strength $\lambda$ on the condition in
Eq.(42), let us consider a case in which Eq.(42) is satisfied at $\lambda=0$
with $w^{\epsilon}_{\rm max}=0$. For example, one has such a case, if the
initial shell $\Gamma^{\mathcal{E}}_{0}$ is sufficiently narrow and the value
of $\Delta_{S}$ is sufficiently small. With increase of $\lambda$ from $0$,
the value of $\delta e$ increases due to the increase of $w^{\epsilon}_{\rm
max}$. At a small $\lambda$, the width $w^{\epsilon}_{\rm max}$ is still small
and, as a result, Eq.(42) is also satisfied.

When the value of $\lambda$ increases beyond some regime, usually, it is
possible for $\Delta h$ to become sufficiently large such that Eq.(42)
gradually becomes invalid. Note that the width $w^{\epsilon}_{\rm max}$ has no
upper bound, because it should increase (approximately) linearly with
$\lambda$ when the interaction Hamiltonian dominates in the total Hamiltonian.
To be quantitative, related to breakdown of Eq.(42), one may consider a value
of $\lambda$, indicated as $\lambda_{h}$, at which the value of $|\Delta
h/h_{0}|$ first reaches $\epsilon_{h}$ when $\lambda$ increases from $0$.
Making use of Eqs.(59) and (36), from Eq.(42) one gets that

$\displaystyle\delta e_{0}+2\Delta_{S}+w^{\epsilon}_{\rm
max}(\lambda_{h})\approx
2N\epsilon_{h}\left|\frac{h_{0}}{\widetilde{h}^{\prime}\left(\frac{e_{0}}{N}\right)}\right|.$
(60)

Two properties are seen from Eq.(60): (1) Since the width $w^{\epsilon}_{\rm
max}$ usually increases with increasing $\lambda$, for systems with $\mu<1$,
the value of $\lambda_{h}$ may increase with increasing $N$; and, (2)
$\lambda_{h}$ should increase with decreasing $\Delta_{S}$, if other
parameters are fixed.

##  IV Further discussions

In this section, we discuss two situations, in which some modified versions of
the RDM relations given in the main result still hold when some restrictions
used above are loosened. In Sec.IV.1, we derive RDM relations in the weak
coupling limit, without the restriction of Eq.(42). In Sec.IV.2, we show that
the main result may be generalized to a generic local interaction Hamiltonian.

###  IV.1 Offdiagonal elements at very weak couplings

In this section, in the weak coupling limit of $\lambda$, without using the
condition in Eq.(42), we derive an expression for offdiagonal elements of the
averaged RDM of nondegenerate levels, by employing a first-order perturbation
treatment. In this limit, diagonal elements of RDM keep approximately
constants, directly given by the initial condition:

$\displaystyle\overline{\rho}^{S}_{\alpha\alpha}\simeq|c_{0\alpha}|^{2}\qquad\text{for
nondegenerate levels $\alpha$}.$ (61)

To be specific, below, we consider two arbitrary nondegenerate levels of the
central system $S$, indicated by $\alpha$ and $\beta$ with $e_{\beta}^{S}\neq
e_{\alpha}^{S}$. The zeroth-order branches, denoted by $|\mathcal{E}^{\rm
0th}_{\alpha}(t)\rangle$, are computed by the Schrödinger evolution of the
initial state $|\Psi(0)\rangle$ under the uncoupled Hamiltonian $H^{0}$.
Noting Eq.(11), one directly gets that

$\displaystyle|\mathcal{E}^{\rm 0th}_{\alpha}(t)\rangle$
$\displaystyle=\langle\alpha|e^{-iH^{0}t}|\Psi(0)\rangle$
$\displaystyle=c_{0\alpha}e^{-ie^{S}_{\alpha}t}\sum_{e_{j}\in\Gamma^{\cal
E}_{0}}e^{-ie_{j}t}c_{0j}|j\rangle.$ (62)

Substituting Eq.(62) into Eq.(21) and noting that $e_{\beta}^{S}\neq
e_{\alpha}^{S}$, one sees that $\overline{\rho}_{\alpha\beta}^{S}$ has a
vanishing zeroth-order term.

The zeroth-order term of $\overline{F}_{\alpha\beta}$, indicated as
$\overline{F}_{\alpha\beta}^{\rm 0th}$, is computed by substituting Eq.(62)
into Eq.(24) and taking the long-time average. Noting that the chaotic
environment $\mathcal{E}$ has a nondegenerate spectrum, direct computation
gives that 777For an environment that possesses a degenerate spectrum, one may
divide the set of those labels $i$, for which
$e_{i}\in\Gamma^{\mathcal{E}}_{0}$, into subsets according to the degeneracy.
We denote the subsets by ${\mathcal{D}}_{q}$ with a label $q$, such that
$e_{i}=e_{j}$ for all $i,j\in{\mathcal{D}}_{q}$. Then, it is easy to find that
$\displaystyle
h_{1}=\sum_{q}\sum_{i,j\in{\mathcal{D}}_{q}}c_{0i}^{*}H^{I\mathcal{E}}_{ij}c_{0j}.$
(63)

$\displaystyle\overline{F}_{\alpha\beta}^{\rm
0th}=|c_{0\alpha}|^{2}h_{1}\delta_{\alpha\beta},$ (64)

where

$\displaystyle
h_{1}=\sum_{e_{i}\in\Gamma^{\mathcal{E}}_{0}}|c_{0i}|^{2}H^{I\mathcal{E}}_{ii}.$
(65)

Now, we compute the first-order term of $\overline{\rho}_{\alpha\beta}^{S}$.
For this purpose, let us rewrite Eq.(29) as follows,

$\displaystyle\overline{\rho}_{\alpha\beta}^{S}=\frac{\lambda}{e_{\beta}^{S}-e_{\alpha}^{S}}\sum_{\gamma=1}^{m}\Big{[}H_{\alpha\gamma}^{IS}\overline{F}_{\beta\gamma}-H^{IS}_{\gamma\beta}\overline{F}_{\gamma\alpha}\Big{]}.$
(66)

Substituting the above-obtained zeroth-order terms
$\overline{F}_{\alpha\beta}^{\rm 0th}$ into the rhs of Eq.(66), one gets the
following expression of $\overline{\rho}_{\alpha\beta}^{S}$ up to the first-
order term:

$\displaystyle\overline{\rho}_{\alpha\beta}^{S}\simeq\frac{\lambda
H_{\alpha\beta}^{IS}h_{1}}{e_{\beta}^{S}-e_{\alpha}^{S}}\left(|c_{0\beta}|^{2}-|c_{0\alpha}|^{2}\right).$
(67)

Finally, we compare two results obtained above, Eq.(67) and Eq.(52), the
latter of which is a TLS case of the main result in Eq.(51). The two results
were gotten under different conditions: Eq.(67) was derived _merely_ under the
condition of very weak coupling, while, Eq.(51) was derived under a condition
that includes three requirements : ETH ansatz in Eq.(10), Eq.(42), and
largeness of $N$. We would remark that the above two conditions are sufficient
conditions for the corresponding results, but not necessary conditions. For
example, it is possible for Eq.(52) to hold in some cases, even when Eqs.(10)
and (42) are not fulfilled. In addition, none of the two conditions includes
the other.

To show consistency of the above two results, let us consider a case in which
both conditions are satisfied. In fact, under Eqs.(10) and (42), it is easy to
see that $h_{1}$ in Eq.(65) satisfies that $h_{1}\simeq h_{0}$. Then, in the
weak coupling limit with Eq.(61), Eq.(67) is written as

$\displaystyle\overline{\rho}^{S}_{\alpha\beta}\simeq\lambda{\eta_{r}h_{0}}(|c_{0\beta}|^{2}-|c_{0\alpha}|^{2}).$
(68)

Clearly, Eq.(68) gives the same prediction as Eq.(52) in this case.

###  IV.2 A generic interaction

In this section, we give a brief discussion for a generic local interaction
Hamiltonian $H^{I}$, which is written as a sum of direct-product terms.
Suppose that there are $M_{\rm LIT}$ such terms, with the subscript “LIT”
standing for “local interaction terms”. Then, $H^{I}$ is written as

$H^{I}=\sum_{\nu=1}^{M_{\rm LIT}}\lambda_{\nu}H^{IS,\nu}\otimes
H^{I\mathcal{E},\nu},$ (69)

where $\lambda_{\nu}$ are parameters and $H^{I\mathcal{E},\nu}$ are local
operators of the environment. The operators $H^{I\mathcal{E},\nu}$ are assumed
to satisfy the ETH ansatz, with functions $h^{\nu}(e)$, respectively. For such
a generic $H^{I}$, the operator $F(t)$ in Eq.(25) is written as

$\displaystyle F(t)=\sum_{\nu=1}^{M_{\rm
LIT}}F_{\alpha\beta}^{\nu}(t)|\alpha\rangle\langle\beta|,$ (70)

where

$F_{\alpha\beta}^{\nu}(t)=\langle\mathcal{E}_{\alpha}(t)|H^{I\mathcal{E},\nu}|\mathcal{E}_{\beta}(t)\rangle.$
(71)

Following arguments similar to those given in Sec.III, with appropriate
generalizations, one may study the long-time average of this generic operator
$F(t)$ and get similar results. More exactly, the main generalization is that
Eq.(42) is now written as

$\displaystyle\left|\frac{1}{h_{0}^{\nu}}\Delta
h^{\nu}\right|\leq\frac{\epsilon_{h}}{M_{\rm LIT}}\quad\text{with
$\epsilon_{h}\ll 1$}\quad(\forall\nu),$ (72)

where $h_{0}^{\nu}=h^{\nu}(e_{0})$ and

$\displaystyle\Delta
h^{\nu}=\max_{e_{i}\in{\Gamma}^{\mathcal{E}}}|{h^{\nu}(e_{i})-h_{0}^{\nu}}|.$
(73)

The final result is that, at a sufficiently large $N$,

$\displaystyle[\widetilde{H}^{S},\overline{\rho}^{S}]\simeq 0,$ (74)

where

$\displaystyle\widetilde{H}^{S}=H^{S}+\sum_{\nu=1}^{M_{\rm
LIT}}\lambda_{\nu}h_{0}^{\nu}H^{IS,\nu}.$ (75)

##  V Numerical tests

In this section, we present numerical simulations that have been performed for
checking analytical predictions given above. Specifically, we discuss the
employed model and analytical predictions in Sec.V.1, and discuss numerical
simulations in Sec.V.2.

###  V.1 The model

In numerical simulations, we employ a TLS as the central system $S$ and one
defect Ising chain as the environment $\mathcal{E}$. The TLS has a self-
Hamiltonian written as

$\displaystyle H^{S}=q_{s}S^{z},$ (76)

where $q_{s}$ is a parameter and $S^{z}$ indicates the $z$-component Pauli
matrix divided by $2$.

The defect Ising chain is composed of a number $N$ of $\frac{1}{2}$-spins
lying in an inhomogeneous transverse field, whose Hamiltonian is written as

$H^{\mathcal{E}}=B_{x}\sum_{l=1}^{N}S_{l}^{x}+d_{1}S_{1}^{z}+d_{5}S_{5}^{z}+J_{z}\sum_{l=1}^{N}S_{l}^{z}S^{z}_{l+1},$
(77)

where $S^{x}_{l}$ and $S^{z}_{l}$ indicate Pauli matrices divided by $2$ at
the $l$-th site. Here, $B_{x}$, $J_{z}$, $d_{1}$, and $d_{5}$ are parameters,
which are adjusted such that the defect Ising chain is a quantum chaotic
system. That is, for levels not close to edges of the energy spectrum, the
nearest-level-spacing distribution $P(s)$ is close to the Wigner-Dyson
distribution $P_{W}(s)=\frac{\pi}{2}s\exp(-\frac{\pi}{4}s^{2})$, the latter of
which is almost identical to the prediction of RMT [44, 45, 46]. Exact values
of the parameters used are $B_{x}=0.9,J_{z}=1.0$, $d_{1}=1.11$, and
$d_{5}=0.6$; and $N$ is between $10$ and $13$. In our numerical computation of
EFs, the periodic boundary condition was implied and the so-called Krylov-
space method was used.

The TLS is coupled to the $k$-th spin of the defect Ising chain. We have
studied two specific forms of the local interaction Hamiltonian, indicated as
$H^{I}_{(1)}$ and $H^{I}_{(2)}$,

$\displaystyle H^{I}_{(1)}=\lambda S^{x}\otimes S_{k}^{x},$ (78a)
$\displaystyle H^{I}_{(2)}=\lambda(S^{x}+S^{z})\otimes S_{k}^{x}.$ (78b)

Their difference lies in that the TLS part of $H^{I}_{(1)}$ has no overlap
with $H^{S}$ in Eq.(76), while, $H^{I}_{(2)}$ has some. According to Eqs.(4)
and (53), one finds that $\Delta_{S}=q_{s}$, $\eta_{d}=0$ and
$\eta_{r}=\frac{1}{2q_{s}}$ for $H^{I}_{(1)}$, and $\eta_{d}=-\frac{1}{q_{s}}$
and $\eta_{r}=\frac{1}{2q_{s}}$ for $H^{I}_{(2)}$. Numerically, we have
checked that the ETH ansatz is applicable to local operators in the defect
Ising chain (see Appendix C).

Figure 2: Values of $|\overline{\rho}^{S}_{12}|$ (triangles) vs the coupling
strength $\lambda$ in the logarithm scale, under two interaction Hamiltonians
of $H^{I}_{(1)}$ and $H^{I}_{(2)}$. Left panels: $q_{s}=0.05$, and right
panels: $q_{s}=0.3$. The solid lines (red) represent predications of the main
result Eq.(52), and the dashed-dotted lines (blue) show predications of
Eq.(67) for very weak couplings. The vertical dot lines (black) indicate
positions of $\lambda_{c}$, which were computed by Eq.(79) with
$\epsilon_{c}=0.1$. Parameters: $(c_{01},c_{02})=(0.51,0.86)$, $N=13$,
$e_{0}=-1.2$, $\delta e_{0}=0.1$, and $k=7$.

Below, we discuss predictions for properties of the long-time-averaged RDM
element $\overline{\rho}^{S}_{12}$ of the TLS, which are given by analytical
results of previous sections. We discuss in the increasing order of the
interaction strength $\lambda$.

(1) Regime of very small $\lambda$ (weak coupling limit).
As discussed in Sec.IV.1, $\overline{\rho}^{S}_{12}$ should satisfy Eq.(67) at
very small $\lambda$. Since the ETH ansatz is applicable to the defect Ising
chain, when Eq.(42) is satisfied, this prediction coincides with Eq.(52),
which is the TLS case of the main result in Eq.(51).

(2) Regime of small but not very small $\lambda$.
(a) Eq.(42) being valid at $\lambda=0$.
In this case, Eq.(42) is also valid at small $\lambda$. As a result,
$\overline{\rho}^{S}_{12}$ should satisfy the main prediction Eq.(52) at
sufficiently large $N$.
(b) Eq.(42) being invalid at $\lambda=0$.
In this case, there is no definite analytical prediction for
$\overline{\rho}^{S}_{12}$ beyond the weak coupling limit.

(3) Regime of $\lambda$ below $\lambda_{h}$ with Eq.(42) valid.
As discussed in Sec.III.4.2, Eq.(52) is applicable for $\lambda$ below
$\lambda_{h}$. The value of $\lambda_{h}$, which satisfies Eq.(60), is
expected to increase with increasing $N$ if $\mu<1$, while, increase with
decreasing $\Delta_{S}$.

As discussed previously, Eq.(42) belongs to a sufficient, but not necessary,
condition for validity of Eq.(52). This implies that Eq.(52) might be useful
even beyond $\lambda_{h}$. To directly study validity of Eq.(52), one may
compute the value of $\lambda$, indicated by $\lambda_{c}$, at which the
relative error first reaches some small parameter indicated by $\epsilon_{c}$
when $\lambda$ increases from $0$,

$\displaystyle\left|\frac{\overline{\rho}^{S}_{12}-\overline{\rho}^{S}_{12,{\rm
th}}}{\overline{\rho}^{S}_{12}}\right|_{\lambda=\lambda_{c}}=\epsilon_{c},$
(79)

where $\overline{\rho}^{S}_{12}$ indicates the exact value of the RDM element
and $\overline{\rho}^{S}_{12,{\rm th}}$ is for the prediction of Eq.(52).

We have no definite analytical prediction for behaviors of $\lambda_{c}$. It
seems reasonable to expect that, at least in some cases, $\lambda_{c}$ may
show some behavior qualitatively similar to that of $\lambda_{h}$ as discussed
above in predication (3).

Figure 3:  (a) Variation of $w_{\rm max}^{\epsilon}$ vs $\lambda$, with
$\lambda$ in the logarithm scale, for $q_{s}=0.05$, $N\in[10,13]$, and
$\epsilon=0.05$. The interaction Hamiltonian is $H_{(1)}^{I}$. Inset: both
axes in the logarithm scale, showing that $w_{\rm max}^{\epsilon}$ has
approximately a $\lambda^{2}$-behavior in the middle regime of $\lambda$. (b)
Similar to (a), but for $q_{s}=0.3$. (c) Variation of $|\Delta h/h_{0}|$ vs
$\lambda$ for $q_{s}=0.05$. Inset: $|\Delta h/h_{0}|$ vs $N$ at $\lambda=0.1$,
showing a $1/N$ behavior, as predicted by Eq.(59). (d) Similar to (c), but for
$q_{s}=0.3$, with some difference in the scale of the vertical axis. Other
parameters are the same as in Fig.2, except that $e_{0}$ is determined by
keeping $e_{0}/N$ constant.

###  V.2 Numerical simulations

We have numerically checked the above predictions for various values of the
parameters concerned. The environmental initial state was taken as a typical
state within an energy shell $\Gamma^{\cal E}_{0}$, which is given by
$e_{0}=-1.2$ and $\delta e_{0}=0.1$.

Two values of $q_{s}$ has been studied, namely, $q_{s}=0.3$ and $0.05$. For
$q_{s}=0.3$, we found that $|\Delta h/h_{0}|\simeq 0.6$ at $\lambda=0$ and
$N=13$, implying invalidity of Eq.(42). With $q_{s}$ changed to $q_{s}=0.05$,
we found $|\Delta h/h_{0}|\simeq 0.1$, implying validity of Eq.(42). In both
cases, $h_{1}\simeq h_{0}$.

Variations of $|\overline{\rho}^{S}_{12}|$ versus the interaction strength
$\lambda$ are shown in Fig.2, for the above-mentioned two values of $q_{s}$
and for the two interaction Hamiltonians in Eq.(78). One sees that there is in
fact no qualitative difference between results for the two interaction
Hamiltonians. In the computation of the rhs of Eq.(52), exact values of
$\overline{\rho}^{S}_{11}$ and $\overline{\rho}^{S}_{22}$ were used. In
agreement with prediction, both the main result of Eq.(52) (solid lines) and
the weak-coupling prediction of Eq.(67) (dashed-dotted lines) work well at
very small $\lambda$, more exactly, at $\lambda$ around $0.001$ and smaller.
Consistently, the mean nearest-level spacing of the total system was found
about $7.3\times 10^{-4}$ in the considered energy region at $N=13$.

With $\lambda$ increased above $0.001$, as expected, the weak-coupling
predictions (dashed-dotted lines in blue) gradually deviate from the exact
values of $|\overline{\rho}^{S}_{12}|$ (triangles). Meanwhile, consistent with
the prediction of (2)(a), for $q_{s}=0.05$ with Eq.(42) valid, predictions of
Eq.(52) (solid lines in red) remain close to the triangles, up to $\lambda\sim
0.1$. It is of interested to note that, even in the case of $q_{s}=0.3$ with
Eq.(42) unsatisfied, predictions of Eq.(52) remains valid up to $\lambda\simeq
0.01$.

To get further understanding for the above-discussed behaviors of
$|\overline{\rho}^{S}_{12}|$, we have studied variation of the maximum width
$w^{\epsilon}_{max}$, which is responsible to the $\lambda$-dependence of the
width $\delta e(=\delta e_{0}+2q_{s}+w_{max}^{\epsilon})$ of
$\Gamma^{\mathcal{E}}$, versus $\lambda$, as well as variation of $|\Delta
h/h_{0}|$ (Fig.3). It is seen that, at $q_{s}=0.05$ and $N=13$, the value of
$w^{\epsilon}_{\rm max}$ keeps small for small $\lambda$ and begin to increase
fast around $\lambda=0.1$; and, consistently, $|\Delta h/h_{0}|$ (triangles
down) behaves in a similar way. Similar behaviors are seen at $q_{s}=0.3$ and
$N=13$, except that $|\Delta h/h_{0}|$ is already large at $\lambda=0$.

We have also studied impact of the particle number $N$. As seen in Fig.3, the
width $w^{\epsilon}_{\rm max}$ is almost independent of $N$ for $N$ from $10$
to $13$, which implies a negligible value of $\mu$, i.e., $\mu\approx 0$.
Meanwhile, the value of $|\Delta h/h_{0}|$ decreases with increasing $N$, in
agreement with a prediction of Eq.(59) that $\Delta h$ may scale as $1/N$ at a
fixed value of $e_{0}/N$ [insets of Fig.3(c) and (d)]. Moreover, with
$\mu\approx 0$, according to prediction (3), $\lambda_{h}$ may increase with
increasing $N$; in other words, Eq.(52) may work better at larger $N$, which
is seen by comparing Fig.4 and Fig.2.

Finally, we discuss numerically obtained values of $\lambda_{h}$ for validity
of Eq.(42) and values of $\lambda_{c}$ for practical use of Eq.(52). The
values of $\lambda_{h}$ may be directly gotten from Fig.3. Taking
$\epsilon_{h}=0.1$, we found that Eq.(42) is valid for no value of $\lambda$
at $q_{s}=0.3$ and $N=10,11,12,13$, which is indicated as $\lambda_{h}={\rm
none}$ in Table 1; and, similarly, for $q_{s}=0.05$ and $N=10$. Thus,
$\lambda_{h}$ has definite values only at $q_{s}=0.05$ and $N=11,12,13$. If
the restriction of $\epsilon_{h}=0.1$ is loosed a little, i.e., if taking
$\epsilon_{h}$ larger than $0.1$ but still small (e.g., $0.15$), $\lambda_{h}$
may have definite values in more cases, which is clear from Fig.3(c). In all
the cases in which $\lambda_{h}$ has definite values for Eq.(42), we found
that larger value of $\lambda_{h}$ corresponds to larger value of $N$,
meanwhile, larger value of $\lambda_{h}$ corresponds to smaller value of
$q_{s}$, in agreement with prediction (3).

Values of $\lambda_{c}$ were computed by making use of Eq.(79) with
$\epsilon_{c}=0.1$. At $q_{s}=0.05$ with Eq.(42) valid, as seen in Table 1,
$\lambda_{c}$ increases with increasing $N$ and is close to $\lambda_{h}$ for
$N=11,12,13$. But, at $q_{s}=0.3$ with Eq.(42) invalid, $\lambda_{c}$ shows a
quite complicated behavior; more exactly, it does not increase monotonically
with $N$ and is unexpectedly large at $N=10$.

| $N=10$ | $N=11$ | $N=12$ | $N=13$
---|---|---|---|---
$\lambda_{h}(q_{s}=0.05)$ | None | 0.04 | 0.07 | 0.1
$\lambda_{h}(q_{s}=0.3)$ | None | None | None | None
$\lambda_{c}(q_{s}=0.05)$ | 0.025 | 0.04 | 0.07 | 0.1
$\lambda_{c}(q_{s}=0.3)$ | 0.1 | 0.015 | 0.025 | 0.015
Table 1: Values of $\lambda_{c}$ and $\lambda_{h}$, obtained with
$\epsilon_{h}=\epsilon_{c}=0.1$, the interaction Hamiltonian is $H_{(1)}^{I}$.
Figure 4: Similar to Fig.2 (a), but, for $N=10$ and $12$.

##  VI Conclusions and discussions

In this paper, the long-time averaged RDM has been studied for a generic small
central system $S$ with $m$ levels, which is locally coupled to a large many-
body chaotic environment $\mathcal{E}$, with the total system undergoing a
Schrödinger evolution. Beside largeness of the particle number of
$\mathcal{E}$, the only restriction is that the environmental part of the
interaction Hamiltonian satisfies the ETH ansatz, with the diagonal term in
the ansatz [namely the function $h(e)$] approximately a constant within the
energy region of relevance. For such a total system, on the eigenbasis of the
central system, $\frac{1}{2}(m-1)(m+2)$ approximate relations have be derived
among elements of its steady states (if existing).

The above-discussed relations imply that the steady RDM should be commutable
with a renormalized Hamiltonian $\widetilde{H}^{S}$ of the central system,
which includes certain averaged impact of the system-environment interaction.
As a consequence, decoherence happens on the eigenbasis of the renormalized
Hamiltonian, even under a system-environment interaction that is dissipative
for the original Hamiltonian $H^{S}$, and leads to a PBS given by the
eigenbasis of $\widetilde{H}^{S}$. This enriches analytical knowledge about
PBS for systems under nonweak and dissipative system-environment interactions,
which had been previously observed numerically in some specific models (see,
e.g., Ref.[47]). 888We would note a difference between the type of models
studied in this paper and the spin-boson models used in Refs.[17, 18, 19, 20,
21]. That is, in a spin-boson model, the environment (the bosons) is not a
quantum chaotic system and the ETH ansatz is usually inapplicable. Due to this
difference, even if a PBS may exist in a spin-boson model, the mechanism
should be quite different from that discussed in this paper.  Moreover,
results of this paper give an explicit way of constructing renormalized
Hamiltonian for PBS.

In fact, renormalized Hamiltonian is also used in a standard master-equation
approach to RDM. There, at an initial stage before the derivation begins, the
self-Hamiltonian of the central system is taken as certain renormalized
Hamiltonian, which we indicate as $H^{S}_{\rm mas}$ with “mas” standing for
“master equation”, given by $H^{S}_{\rm mas}=H^{S}+\lambda
H^{IS}\text{tr}(H^{I\mathcal{E}}\rho^{\mathcal{E}}_{\rm th})$, where
$\rho^{\mathcal{E}}_{\rm th}$ denotes a thermal state of the environment.
Under the ETH ansatz and the condition in Eq.(42), $H^{S}_{\rm mas}$ has
almost the same expression as $\widetilde{H}^{S}$ in Eq.(50), if the state
$\rho^{\mathcal{E}}_{\rm th}$ lies effectively within the energy shell
$\Gamma^{\mathcal{E}}$.

However, there is a big difference between the physical meanings of
$H^{S}_{\rm mas}$ and $\widetilde{H}^{S}$. In fact, in our approach, the
operator $\widetilde{H}^{S}$ is derived by faithfully taking the long-time
average over the overall Schrödinger evolution; and it indicates the existence
of a PBS, if the RDM may approach a steady state. While, in the master-
equation approach, the operator $H^{S}_{\rm mas}$ is mainly employed for the
sake of convenience in derivation, though with deep physical intuition lying
behind it. Only after a certain type of analytical solution to a derived
master equation is found, which is usually a hard task except in some special
models, could it become clear whether $H^{S}_{\rm mas}$ may indeed be of
relevance to a PBS. Moreover, as an approach based a perturbative treatment,
validity of the master-equation approach at long times is a subtle issue.

Finally, we would mention that, beside the field of decoherence, results of
this paper may also be useful in other fields in which properties of steady
states of small and open quantum systems are of relevance, such as quantum
thermodynamics [48, 49, 36, 50].

###### Acknowledgements.

This work was partially supported by the Natural Science Foundation of China
under Grant No. 11535011, 11775210, and 12175222. JW are supported by the
Deutsche Forschungsgemeinschaft (DFG) within the Research Unit FOR 2692 under
Grant No. 397107022 (GE 1657/3-2).

##  Appendix A Derivation of Eq.(22)

In this appendix, we derive Eq.(22). Using Eq.(21), the time evolution of the
RDM is written as

$i\hbar\frac{d\rho_{\alpha\beta}^{S}(t)}{dt}=i\hbar\frac{d}{dt}\langle\mathcal{E}_{\beta}(t)|\mathcal{E}_{\alpha}(t)\rangle=A_{1}+A_{2},$
(80)

where

$\displaystyle A_{1}$
$\displaystyle=i\hbar\left(\frac{d}{dt}\langle\mathcal{E}_{\beta}(t)|\right)|\mathcal{E}_{\alpha}(t)\rangle,$
(81a) $\displaystyle A_{2}$
$\displaystyle=i\hbar\langle\mathcal{E}_{\beta}(t)|\left(\frac{d}{dt}|\mathcal{E}_{\alpha}(t)\rangle\right).$
(81b)

Making use of Eq.(19), one finds that

$\displaystyle
A_{1}=-\sum_{\gamma}\langle\mathcal{E}_{\gamma}(t)|H_{\gamma\beta}|\mathcal{E}_{\alpha}(t)\rangle.$
(82)

From Eqs.(7) and (20), one gets that

$\displaystyle
H_{\alpha\beta}=\langle\alpha|H|\beta\rangle=e^{S}_{\alpha}\delta_{\alpha\beta}+\lambda
H^{IS}_{\alpha\beta}H^{I\mathcal{E}}+H^{\mathcal{E}}\delta_{\alpha\beta}.$
(83)

Then, we write Eq.(82) as,

$\displaystyle A_{1}=$ $\displaystyle-
e_{\beta}^{S}\langle\mathcal{E}_{\beta}(t)|\mathcal{E}_{\alpha}(t)\rangle-\lambda\sum_{\gamma}H^{IS}_{\gamma\beta}\langle\mathcal{E}_{\gamma}(t)|H^{I\mathcal{E}}|\mathcal{E}_{\alpha}(t)\rangle$
$\displaystyle-\langle\mathcal{E}_{\beta}(t)|H^{\mathcal{E}}|\mathcal{E}_{\alpha}(t)\rangle.$
(84)

Noting Eqs.(21) and (24), the above equality gives that

$\displaystyle
A_{1}=-e_{\beta}^{S}\rho^{S}_{\alpha\beta}(t)-\lambda\sum_{\gamma}H^{IS}_{\gamma\beta}F_{\gamma\alpha}(t)-\langle\mathcal{E}_{\beta}(t)|H^{\mathcal{E}}|\mathcal{E}_{\alpha}(t)\rangle.$
(85)

Similarly, one finds

$A_{2}=e_{\alpha}^{S}\rho_{\alpha\beta}^{S}(t)+\lambda\sum_{\gamma}H^{IS}_{\alpha\gamma}F_{\beta\gamma}(t)+\langle\mathcal{E}_{\beta}(t)|H^{\mathcal{E}}|\mathcal{E}_{\alpha}(t)\rangle.$
(86)

Putting the above results together, one gets Eq.(22).

##  Appendix B Scaling of the fluctuation operator

In this appendix, we show that the main $N$-scaling behavior of the
fluctuation operator $\Delta$,

$\Delta=\sum_{\alpha\beta}\Delta_{\alpha\beta}|\alpha\rangle\langle\beta|,$
(87)

is an exponential decay with increasing $N$. For this purpose, let us compute
the Frobenius norm of $\Delta$,

$\displaystyle||\Delta||_{F}^{2}$
$\displaystyle=\sum_{\alpha\beta}|\Delta_{\alpha\beta}|^{2}.$ (88)

Making use of Eq.(41), direct derivation shows that

$\displaystyle||\Delta||_{F}^{2}$
$\displaystyle=\sum_{\alpha\beta}\left|\sum_{ij}g(e,\omega)e^{-S(e)/2}\overline{f_{\beta
i}^{*}f_{\alpha j}}R_{ij}\right|^{2}$
$\displaystyle=\sum_{\alpha\beta}\sum_{iji^{\prime}j^{\prime}}g(e,\omega)e^{-S(e)/2}\overline{f_{\beta
i}^{*}f_{\alpha j}}R_{ij}$
$\displaystyle\quad\quad\times{g}^{*}(e^{\prime},\omega^{\prime})e^{-S(e^{\prime})/2}\overline{f_{\beta
i^{\prime}}f_{\alpha j^{\prime}}^{*}}\ {R}^{*}_{i^{\prime}j^{\prime}}.$ (89)

Note that $g(e,\omega)={g}^{*}(e,-\omega)$ and $R_{ij}={R}^{*}_{ji}$. To
proceed, let us discuss the statistical average of $||\Delta||_{F}^{2}$, taken
over the random variables $R_{ij}$, which is indicate by
$\langle\cdot\rangle$. This averaging procedure results in that [51, 26]

$\displaystyle\langle
R_{ij}{R}^{*}_{i^{\prime}j^{\prime}}\rangle=\delta_{ii^{\prime}}\delta_{jj^{\prime}}+\delta_{ij^{\prime}}\delta_{i^{\prime}j},$
(90)

and, as a consequence,

$\displaystyle\langle||\Delta||_{F}^{2}\rangle$
$\displaystyle=\sum_{ij}|g(e,\omega)|^{2}e^{-S(e)}$
$\displaystyle\quad\quad\times\left(\sum_{\alpha\beta}\big{|}\overline{f_{\beta
i}^{*}f_{\alpha j}}\big{|}^{2}+\overline{f_{\beta i}^{*}f_{\alpha
j}}\cdot\overline{f_{\beta j}f_{\alpha i}^{*}}\right)$ $\displaystyle\leq
2\sum_{ij}|g(e,\omega)|^{2}e^{-S(e)}\sum_{\alpha\beta}\big{|}\overline{f_{\beta
i}^{*}f_{\alpha j}}\big{|}^{2}$ $\displaystyle\leq
2\max_{ij}(|g(e,\omega)|^{2}e^{-S(e)})\sum_{ij}\sum_{\alpha\beta}\big{|}\overline{f_{\beta
i}^{*}f_{\alpha j}}\big{|}^{2}.$ (91)

To compute $\overline{f_{\beta i}^{*}f_{\alpha j}}$, we make use of the fact
that $f_{\alpha j}=\langle\alpha j|\Psi(t)\rangle$. This gives that

$\displaystyle\overline{f_{\beta i}^{*}f_{\alpha
j}}=\overline{\langle\Psi(t)|\beta i\rangle\langle\alpha j|\Psi(t)\rangle}$
$\displaystyle=\sum_{mn}\langle\Psi(0)|n\rangle\langle n|\beta
i\rangle\overline{e^{i(E_{n}-E_{m})t}}\langle\alpha j|m\rangle\langle
m|\Psi(0)\rangle.$ (92)

Note that the environment $\mathcal{E}$, as a quantum chaotic system, has a
nondegenerate spectrum. Under a generic system-environment interaction, the
spectrum of the total system is nondegenerate, too. Then, one has
$\overline{e^{i(E_{n}-E_{m})t}}=\delta_{mn}$ and, as a result,

$\displaystyle\overline{f_{\beta i}^{*}f_{\alpha
j}}=\sum_{n}\langle\Psi(0)|n\rangle\langle n|\beta i\rangle\langle\alpha
j|n\rangle\langle n|\Psi(0)\rangle.$ (93)

This gives that

$\displaystyle\left|\overline{f_{\beta i}^{*}f_{\alpha j}}\right|^{2}$
$\displaystyle=\sum_{nm}\langle\Psi(0)|n\rangle\langle n|\beta
i\rangle\langle\beta i|m\rangle\langle m|\Psi(0)\rangle$
$\displaystyle\times\langle\Psi(0)|m\rangle\langle m|\alpha
j\rangle\langle\alpha j|n\rangle\langle n|\Psi(0)\rangle.$ (94)

Then, making use of the completeness of the basis of $|\alpha\rangle$ and that
of $|i\rangle$, one gets that

$\displaystyle\sum_{\alpha\beta ij}\left|\overline{f_{\beta i}^{*}f_{\alpha
j}}\right|^{2}=\frac{1}{L_{0}},$ (95)

where $L_{0}$ is the so-called participation function of the initial state
$|\Psi(0)\rangle$, defined by

$\displaystyle
L_{0}=\frac{1}{\sum_{n}\left|\langle\Psi(0)|n\rangle\right|^{4}}.$ (96)

As is known, $L_{0}$ gives a measure to the localization length, i.e., to the
number of those levels $E_{n}$ that are effectively occupied by the state
$|\Psi(0)\rangle$. For a large environment and an initial shell not extremely
narrow, the value of $L_{0}$ is large.

Substituting Eq.(95) into Eq.(91), we get an upper bound to the averaged norm
$\langle||\Delta||_{F}^{2}\rangle$, i.e.,

$\displaystyle\langle||\Delta||_{F}^{2}\rangle$
$\displaystyle\leq\frac{2}{L_{0}}\max_{ij}|g(e,\omega)|^{2}e^{-S(e)}\sim\frac{2N^{2\gamma}}{L_{0}}e^{-S(e)}.$
(97)

Since the averaging procedure does not change the $N$-scaling behavior of the
norm $||\Delta||_{F}^{2}$ and the exponential-decay term $e^{-S(e)/2}$ already
exists in the exact expression of $\Delta_{\alpha\beta}$ in Eq.(41), from
Eq.(97) one sees that the $N$-scaling behavior of fluctuation operator
$\Delta$ should be dominated by the exponential decay $e^{-S(e)/2}$.

##  Appendix C Verification of ETH ansatz

Figure C. 1:  Upper panel: Diagonal matrix elements of two local observables
$S_{k}^{z}$ and $S_{k}^{x}$ ($k=7$) of the defect Ising chain vs the
environmental energy $\epsilon_{i}=e_{i}/N$ for different chain size $N$. In
agreement with the ETH ansatz in Eq.(10), these elements fluctuate around
certain slowly-varying functions of $e$, respectively, and the fluctuations
decrease with the increase of the particle number $N$. Lower panel: Locally
averaged values of the above elements (within windows with a width $0.01$),
showing a feature of approximate size-independence.

Due to the hypothesis feature of the ETH ansatz in Eq.(10), we have checked
its validity in the model employed in this paper. We did this for the two
local operators $S_{k}^{x}$ and $S_{k}^{z}$ at the site $k=7$ in the defect
Ising chain.

Figure C. 2:  (a) Exponential decay of the deviation $\sigma_{d}^{a}$ in
Eq.(99) with the increase of $N$, for fluctuations of the diagonal elements of
$S_{k}^{x}$ (empty squares) and of $S_{k}^{z}$ (solid circles). (b)
Exponential decay of $\sigma^{a}_{nd}$ in Eq.(101) for fluctuations of
offdiagonal elements. The results are in agreement with the prediction of ETH
in Eq.(10).

_Diagonal ETH_.— Let us first discuss predictions of Eq.(10) for diagonal
elements of local observables. Expectation values of the two local
observables,

$\displaystyle(S_{k}^{a})_{ii}=\langle i|S_{k}^{a}|i\rangle\quad\text{with
$a=x,z$,}$ (98)

are plotted in Fig.C. 1. It is seen that, in agreement with ETH, the diagonal
elements fluctuate around certain slowly varying function $h(e)$ and the
fluctuations decrease with increasing the particle number $N$. Note that the
horizontal axis is labeled by $e_{i}/N$. For $a=z$, the values of $h(e)$ are
close to zero, while, for $a=x$, most of $|h(e)|$ are notably larger than
zero.

To study quantitatively the fluctuations of $(S_{k}^{a})_{ii}$, we have
computed the standard deviations $\sigma^{a}_{d}$,

$\sigma_{d}^{a}=\sqrt{\frac{1}{N_{\Gamma_{0}^{\mathcal{E}}}}\sum_{e_{i}\in\Gamma_{0}^{\mathcal{E}}}|(S_{k}^{a})_{ii}-\mu^{a}|^{2}},$
(99)

where

$\mu^{a}=\frac{1}{N_{\Gamma_{0}^{\mathcal{E}}}}\sum_{e_{i}\in\Gamma_{0}^{\mathcal{E}}}(S_{k}^{a})_{ii}.$
(100)

As seen in Fig. C. 2 (a), the fluctuation decays exponentially with the
increase of $N$, as predicted by the term $e^{-S(e)}$ in the second part on
the rhs of Eq.(10). Moreover, in agreement with the prediction of ETH, the
distributions of $[(S_{k}^{a})_{ii}-\mu^{a}]/\sigma^{a}_{d}$ are close to the
Gaussian form [Fig. C. 3 (a) and (b)].

Figure C. 3:  Distributions of fluctuations of the diagonal elements of
$S_{k}^{z}$ (a) and of $S_{k}^{x}$ (b), rescaled by $\sigma_{d}^{z}$ and
$\sigma_{d}^{x}$, respectively. And, distributions of the offdiagonal elements
of $S_{k}^{z}$ (c) and of $S_{k}^{x}$ (d), rescaled by $\sigma_{nd}^{z}$ and
$\sigma_{nd}^{x}$, respectively. The dashed curves represent the Gaussian
distribution with unit variance.  Figure C. 4: Locally averaged values (in the
logarithm scale) of the absolute square of off-diagonal elements of (a)
$S_{k}^{z}$ and (b) $S_{k}^{x}$, within an energy shell centered at $-1.2$ and
with a width $0.2$, vs $\omega=e_{j}-e_{i}$. Local averages were taken within
small windows with width $0.01$.

_Off-diagonal ETH_.— Next, we discuss the offdiagonal elements
$(S_{k}^{a})_{ij}$. In agreement with the prediction of ETH, the probability
distributions of $(S_{k}^{a})_{ij}/\sigma^{a}_{nd}$ have a Gaussian form [Fig.
C. 3 (c) and (d)], where $\sigma^{a}_{nd}$ are the standard deviations for the
offdiagonal elements,

$\sigma_{nd}^{a}=\sqrt{\frac{1}{N_{\Gamma_{0}^{\mathcal{E}}}(N_{\Gamma_{0}^{\mathcal{E}}}-1)}\sum_{i\neq
j\in\Gamma_{0}^{\mathcal{E}}}|(S_{k}^{a})_{ij}|^{2}}.$ (101)

These standard deviations also decay exponentially with the increase of $N$
[Fig. C. 2 (b)].

To get some knowledge about shapes of the function $g(e,\omega)$, which lacks
an analytical expression, numerical simulations have been performed for
locally averaged values of $\overline{|(S_{k}^{z})_{ij}|^{2}}$ and
$\overline{|(S_{k}^{x})_{ij}|^{2}}$ for off-diagonal elements. As seen in
Fig.C. 4, the function shows a size-independent feature, with an exponential-
type decay at large $\omega$.

## References

  * Leggett  _et al._ [1987] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. Fisher, A. Garg, and W. Zwerger, Dynamics of the dissipative two-state system, Reviews of Modern Physics 59, 1 (1987).
  * Breuer  _et al._ [2002] H.-P. Breuer, F. Petruccione, _et al._ ,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * Alicki and Lendi [2007] R. Alicki and K. Lendi,  _Quantum dynamical semigroups and applications_ , Vol. 717 (Springer, 2007).
  * Breuer  _et al._ [2016] H.-P. Breuer, E.-M. Laine, J. Piilo, and B. Vacchini, Colloquium: Non-markovian dynamics in open quantum systems, Reviews of Modern Physics 88, 021002 (2016).
  * De Vega and Alonso [2017] I. De Vega and D. Alonso, Dynamics of non-markovian open quantum systems, Reviews of Modern Physics 89, 015001 (2017).
  * Zurek [1981] W. H. Zurek, Pointer basis of quantum apparatus: Into what mixture does the wave packet collapse?, Physical review D 24, 1516 (1981).
  * Zurek [2003] W. H. Zurek, Decoherence, einselection, and the quantum origins of the classical, Reviews of modern physics 75, 715 (2003).
  * Schlosshauer [2005] M. Schlosshauer, Decoherence, the measurement problem, and interpretations of quantum mechanics, Reviews of Modern physics 76, 1267 (2005).
  * Wiseman and Milburn [2009] H. M. Wiseman and G. J. Milburn,  _Quantum measurement and control_ (Cambridge university press, 2009).
  * Paz and Zurek [1999] J. P. Paz and W. H. Zurek, Quantum limit of decoherence: Environment induced superselection of energy eigenstates, Physical Review Letters 82, 5181 (1999).
  * Joos  _et al._ [2013] E. Joos, H. D. Zeh, C. Kiefer, D. J. Giulini, J. Kupsch, and I.-O. Stamatescu,  _Decoherence and the appearance of a classical world in quantum theory_ (Springer Science & Business Media, 2013).
  * Gorin  _et al._ [2004] T. Gorin, T. Prosen, T. Seligman, and W. Strunz, Connection between decoherence and fidelity decay in echo dynamics, Physical Review A 70, 042105 (2004).
  * Albash and Lidar [2015] T. Albash and D. A. Lidar, Decoherence in adiabatic quantum computation, Physical Review A 91, 062320 (2015).
  * Weiss [2012] U. Weiss,  _Quantum dissipative systems_ , Vol. 13 (World scientific, 2012).
  * Wang  _et al._ [2008] W.-g. Wang, J. Gong, G. Casati, and B. Li, Entanglement-induced decoherence and energy eigenstates, Physical Review A 77, 012108 (2008).
  * He and Wang [2014] L. He and W.-g. Wang, Statistically preferred basis of an open quantum system: Its relation to the eigenbasis of a renormalized self-hamiltonian, Physical Review E 89, 022125 (2014).
  * Lee  _et al._ [2012] C. K. Lee, J. Cao, and J. Gong, Noncanonical statistics of a spin-boson model: Theory and exact monte carlo simulations, Physical Review E 86, 021109 (2012).
  * Addis  _et al._ [2014] C. Addis, G. Brebner, P. Haikka, and S. Maniscalco, Coherence trapping and information backflow in dephasing qubits, Physical Review A 89, 024101 (2014).
  * Roszak  _et al._ [2015] K. Roszak, R. Filip, and T. Novotnỳ, Decoherence control by quantum decoherence itself, Scientific reports 5, 1 (2015).
  * Zhang  _et al._ [2015] Y.-J. Zhang, W. Han, Y.-J. Xia, Y.-M. Yu, and H. Fan, Role of initial system-bath correlation on coherence trapping, Scientific reports 5, 1 (2015).
  * Guarnieri  _et al._ [2018] G. Guarnieri, M. Kolář, and R. Filip, Steady-state coherences by composite system-bath interactions, Physical review letters 121, 070401 (2018).
  * Deutsch [1991] J. M. Deutsch, Quantum statistical mechanics in a closed system, Physical review a 43, 2046 (1991).
  * Srednicki [1994] M. Srednicki, Chaos and quantum thermalization, Physical review e 50, 888 (1994).
  * Srednicki [1999] M. Srednicki, The approach to thermal equilibrium in quantized chaotic systems, Journal of Physics A: Mathematical and General 32, 1163 (1999).
  * Rigol  _et al._ [2008] M. Rigol, V. Dunjko, and M. Olshanii, Thermalization and its mechanism for generic isolated quantum systems, Nature 452, 854 (2008).
  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics, Advances in Physics 65, 239 (2016).
  * Deutsch [2018] J. M. Deutsch, Eigenstate thermalization hypothesis, Reports on Progress in Physics 81, 082001 (2018).
  * Garrison and Grover [2018] J. R. Garrison and T. Grover, Does a single eigenstate encode the full hamiltonian?, Physical Review X 8, 021026 (2018).
  * Wang  _et al._ [2022] J. Wang, M. H. Lamann, J. Richter, R. Steinigeweg, A. Dymarsky, and J. Gemmer, Eigenstate thermalization hypothesis and its deviations from random-matrix theory beyond the thermalization time, Physical Review Letters 128, 180601 (2022).
  * Khatami  _et al._ [2013] E. Khatami, G. Pupillo, M. Srednicki, and M. Rigol, Fluctuation-dissipation theorem in an isolated system of quantum dipolar bosons after a quench, Physical review letters 111, 050403 (2013).
  * Abanin  _et al._ [2015] D. A. Abanin, W. De Roeck, and F. Huveneers, Exponentially slow heating in periodically driven many-body systems, Physical review letters 115, 256803 (2015).
  * Mukerjee  _et al._ [2006] S. Mukerjee, V. Oganesyan, and D. Huse, Statistical theory of transport by strongly interacting lattice fermions, Physical Review B 73, 035113 (2006).
  * Brenes  _et al._ [2020a] M. Brenes, T. LeBlond, J. Goold, and M. Rigol, Eigenstate thermalization in a locally perturbed integrable system, Physical review letters 125, 070605 (2020a).
  * LeBlond  _et al._ [2020] T. LeBlond, D. Sels, A. Polkovnikov, and M. Rigol, Universality in the onset of quantum chaos in many-body systems, arXiv preprint arXiv:2012.07849 (2020).
  * Luitz and Lev [2016] D. J. Luitz and Y. B. Lev, Anomalous thermalization in ergodic systems, Physical review letters 117, 170404 (2016).
  * Wang [2012] W.-g. Wang, Statistical description of small quantum systems beyond the weak-coupling limit, Physical Review E 86, 011115 (2012).
  * [37] W.-g. Wang, Q.-c. Li, M. Yuan, and J. Wang, in preparation.
  * Wang and Wang [2019] J. Wang and W.-g. Wang, Convergent perturbation expansion of energy eigenfunctions on unperturbed basis states in classically-forbidden regions, Journal of Physics A: Mathematical and Theoretical 52, 235204 (2019).
  * Wang  _et al._ [1998] W.-g. Wang, F. Izrailev, and G. Casati, Structure of eigenstates and local spectral density of states: A three-orbital schematic shell model, Physical Review E 57, 323 (1998).
  * Wang [2000] W.-g. Wang, Perturbative and nonperturbative parts of eigenstates and local spectral density of states: The wigner-band random-matrix model, Physical Review E 61, 952 (2000).
  * Wang [2002] W.-g. Wang, Nonperturbative and perturbative parts of energy eigenfunctions: A three-orbital schematic shell model, Physical Review E 65, 036219 (2002).
  * Kim  _et al._ [2014] H. Kim, T. N. Ikeda, and D. A. Huse, Testing whether all eigenstates obey the eigenstate thermalization hypothesis, Physical Review E 90, 052105 (2014).
  * Brenes  _et al._ [2020b] M. Brenes, J. Goold, and M. Rigol, Low-frequency behavior of off-diagonal matrix elements in the integrable xxz chain and in a locally perturbed quantum-chaotic xxz chain, Physical Review B 102, 075127 (2020b).
  * Haake [2013] F. Haake,  _Quantum Signatures of Chaos_ , Vol. 54 (Springer Science & Business Media, 2013).
  * Casati  _et al._ [1980] G. Casati, F. Valz-Gris, and I. Guarnieri, On the connection between quantization of nonintegrable systems and statistical theory of spectra, Lettere al Nuovo Cimento 28, 279 (1980).
  * Bohigas  _et al._ [1984] O. Bohigas, M.-J. Giannoni, and C. Schmit, Characterization of chaotic quantum spectra and universality of level fluctuation laws, Physical review letters 52, 1 (1984).
  * Wang  _et al._ [2012] W.-g. Wang, L. He, and J. Gong, Preferred states of decoherence under intermediate system-environment coupling, Physical review letters 108, 070403 (2012).
  * Binder  _et al._ [2019] F. Binder, L. A. Correa, C. Gogolin, J. Anders, and G. Adesso,  _Thermodynamics in the quantum regime: fundamental aspects and new directions_ , Vol. 195 (Springer, 2019).
  * Goold  _et al._ [2016] J. Goold, M. Huber, A. Riera, L. Del Rio, and P. Skrzypczyk, The role of quantum information in thermodynamics—a topical review, Journal of Physics A: Mathematical and Theoretical 49, 143001 (2016).
  * Wang [2018] W.-g. Wang, Decoherence approach to energy transfer and work done by slowly driven systems, Physical Review E 97, 012128 (2018).
  * Murthy and Srednicki [2019] C. Murthy and M. Srednicki, Bounds on chaos from the eigenstate thermalization hypothesis, Physical review letters 123, 230606 (2019).