# A Compact Formulation for the Multi-UAV Assisted Data Gathering Problem in
Wireless Sensor Networks : A Linear Programming Approach

Chi-Hieu Nguyen  Hanoi University of Science and Technology
Ha Noi, Viet Nam
<EMAIL_ADDRESS>Khanh-Van Nguyen  Hanoi University of Science
and Technology
Ha Noi, Viet Nam
<EMAIL_ADDRESS>
###### Abstract

In this paper, we study the problem of gathering data from large-scale
wireless sensor networks using multiple unmanned air vehicles (UAVs) to gather
data at designated rendezvouses, where the goal is to maximize the network
lifetime. Previous proposals often consider a practical approach where the
problem of determining a data gathering scheme is decomposed into 2 sub-
problems: i) partitioning the networks into clusters for determining the
rendezvouses as these obtained cluster heads; and ii) determining the paths
for a set of a given number of UAVs to come gathering data at these
rendezvouses which have been harvesting data within each local clusters,
respectively. We try to deal with this as a whole optimization problem,
expecting a significant increase in computation complexity which would bring
new challenge in creating practical solutions for large-scale WSNs. We
introduce two alternatives mixed-integer linear programming (MILP)
formulations, namely the 2-index model with $O(n^{2})$ variables and the
3-index model that has $O(n^{3})$ variables, where $n$ denotes the number of
sensor nodes. We show that our best model could solve optimally the problem
instances with up to 50 sensor nodes in less than 30 minutes. Next, we propose
a heuristic idea to reduce the number of variables in implementing the 3-index
model to effectively handle larger-scale networks with size in hundreds. The
experiment results show that our heuristic approach significantly prolongs the
network lifetime compared to existing most efficient proposals.

##  I Introduction

Wireless sensor networks (WSNs) play an important role in the Internet of
Things (IoT) for implementing real-time monitoring or long-term surveillance
purposes. In large-scale WSN, several thousands of sensor nodes powered by
batteries are deployed in a vast region and they coordinate to perform two
main tasks: sensing and communicating. Traditionally, the sensory data must be
periodically sent to a static sink node through multi-hops communication.
However, this approach is suffering from the hot spot problem where the energy
of nodes nearer to the sink depletes faster due to the huge traffic load from
far away nodes [1]. To tackle this issue, employing one or more mobile data
collectors as the mobile sinks to move around and gather sensory data has been
considered as an efficient approach to reduce and also balance the energy
consumption between SNs thus prolong the network lifetime [2].

Recently, the Unmanned Aerial Vehicles (UAVs) technology have seen rapid
growth with application in various fields such as military, emergency rescue,
transportation, surveillance and monitoring. This technology has also emerged
as a flexible and convenient approach to collect the harvested info data in
large-scale WSNs because of their advantages such as low cost, simple
deployment and ability to avoid terrain obstacles compared with other types of
mobile collectors. In additional, the ground-to-air communication links
between sensor node (SN) and UAVs typically exploit the line-of-sight (LoS)
propagation thus the data transmission is more reliable and also provide
higher data rate compared to ground vehicles. However, the critical issue in
this technology is that the UAVs usually have limited endurance due to the
practical physical constraints, and need to be collected for battery swap or
recharging, e.g., the endurance time can be 30 minutes for a typical rotary-
wing UAV [3].

Similar to other kinds of mobile collector, recent research works [4, 5, 6, 7,
8] in UAV-assisted data gathering can be classified into two basic approaches:
rendezvous-less and rendezvous-based. The former approach requires the UAVs to
visit every SN neighborhood at least once so that each SN can send its data
directly to a UAV via single-hop transmission. Although the rendezvous-less
approach improves the communication reliability and energy conservation of
SNs, the length of UAV trajectory becomes excessively long thus this approach
is often infeasible for large-scale WSNs. This challenging issue can be
resolved by using the rendezvous-based approach, where the UAVs only visit a
subset of predefined nodes called rendezvous nodes to gather data from their
nearby SNs.

Several works [9, 10, 7, 11] on the rendezvous-based data gathering have been
proposed and most of them follow a common idea of decomposing the main problem
into two sub-problems which will be solved one by one: i) designing a
clustering scheme wherein the cluster heads (CHs) take the responsibility as
rendezvous nodes to gather data from each cluster; ii) determining UAV
trajectories for collecting the data at the CHs and transporting to the base
sink. Although these sub-problems are easier to deal with, this method can
affect the solution quality.

In this paper, we address at the same time the optimization problem of both
the network clustering and UAV trajectory designing for a data collection
period (or a round in other words) in WSN. Our objective is to maximize the
network lifetime, which can be employed in each round by using an alternative
objective function of SN energy. To this end, we first give a formulation of
the problem under a mixed-integer linear programming model. Although there
exist previous works in the literature that study the UAV-assisted data
gathering problem using linear programming approach, however, they are not
general enough that the proposed model considers only the scenario of a single
UAV as in [12, 13, 14, 15] or follow the rendezvous-less approach as in [16].
To the best of our knowledge, this work is the first attempting to propose a
MILP formulation for the rendezvous-based data gathering problem with multiple
UAVs in WSN. Our contributions are as follows:

  * •

We introduce two formulations of the problem under a mixed-integer linear
programming model, namely the 2-index and 3-index models. Then we show that
the 3-index model can bring significant savings in computational effort
compared to the 2-index model, even though it has a higher number of variables
and constraints.

  * •

We propose a heuristic model that reduces the number of variables in the
3-index model to solve the problem instance representing a practical larger-
scale network with size in hundreds.

  * •

We conduct experiments to demonstrate our heuristic model’s effectiveness in
terms of extending network lifetime compared to the existing proposal.

The rest of this paper is organized as follows. Section II describes the
system model. Section III presents two proposed MILP formulations for the
problem. Section IV discusses our numerical and experimental results. Section
V concludes the paper.

##  II System model

We consider a sensor network consisting of a sink node $s_{0}$ and a set of
$n$ sensor nodes (SNs) $\mathcal{S}=\\{s_{1},...,s_{n}\\}$. Assume that the
data of all sensors must be collected at the sink periodically with a
completion time required of $T_{max}$ seconds, which is determined from the
application requirements. Each period of data collection is called a round and
$T_{max}$ is called the round length. There are $m$ rotate wing UAVs deployed
to collect data from sensors, all fly at a fixed altitude $h$. In each round,
the network is partitioned into clusters and each SN transmits its data to an
associated CH in a multi-hop fashion. Then each UAV departs from the sink,
sequentially visits some specified CHs to collect data and finally return to
the sink. Assume that the residual energy of SN $s_{i}$ after the $r^{th}$
round is $E_{i}^{r}$ the amount of data to collect at $s_{i}$ is $b_{i}$. The
problem asks for a data collection scheme in the $(r+1)^{th}$ round including
$k\geq 1$ cluster trees and $m$ UAV trajectories so that the minimum SN energy
after $(r+1)^{th}$ round is maximized, while satisfying the following:

  * •

Each cluster tree is rooted at the CH and spans all cluster nodes.

  * •

Each SN is included in exactly one cluster tree.

  * •

Each UAV trajectory is a ring that contains the sink.

  * •

Each CH is included in exactly one UAV trajectory.

  * •

Each UAV trajectory has length less than a given value $L_{max}$

We utilize the same radio model as discussed in [17]. In this model, the free
space ($d^{2}$ power loss) model is used if the distance between the
transmitter and the receiver is less than a $d_{0}$ threshold, otherwise the
model of multi-path ($d^{4}$ power loss) is used. Thus, the energy spent to
transmit a $l$-bits packet over distance $d$ is:

$E_{Tx}(l,d)=\begin{cases}lE_{elec}+l\epsilon_{fs}d^{2},&d<d_{0}\\\
lE_{elec}+l\epsilon_{mp}d^{4},&d\geq d_{0}\end{cases}$ (1)

The energy spent by the radio to receive $l$-bits of data is:

$E_{Rx}(l)=lE_{elec}$ (2)

For more convenient, let $E_{Rx}=E_{elec}$ be the per-bit receiving energy of
a SN and $E_{Tx}^{ij}=E_{elec}+\epsilon_{fs}d_{ij}^{2}$ be the per-bit
transmitting energy from node $s_{i}$ to $s_{j}$, respectively. Thus, if node
$s_{i}$ has to relay $B_{i}$ bits of data from the descendants in its cluster
tree in a round, then the total energy consumption of $s_{i}$ after the round
is calculated as:

$\begin{split}E_{i}^{c}&=E_{Rx}(B_{i})+E_{Tx}(B_{i}+b_{i})\\\
&=(E_{Rx}+E_{Tx}^{ij})B_{i}+E_{Tx}^{ij}b_{i}\end{split}$ (3)

The remaining energy of $s_{i}$ after the $(r+1)^{th}$ round is

$E_{i}^{r+1}=E_{i}^{r}-E_{i}^{c}$ (4)

The multi-UAVs data collection problem asks for $m$ data forwarding trees
corresponding to $m$ clusters, each of which roots at a CH node and $k$ UAV
trajectories which start and end at the sink. The objective is to maximize the
network lifetime, which is defined as the time span from deployment to the
instant when the first SN run out of energy.

##  III MILP Formulations

###  III-A Objective function and basic constraints

Both formulations proposed in this section used the same objective function
and shared some basic variables and constraints. We will first present these
fundamentals in this subsection.

In both formulations, the cluster configuration (i.e the network spanning
forest) and all UAV trajectories will be defined by two types of binary
variables $x_{ij}$ and $r_{ij}$. In particular, the variable $x_{ij}$ takes
value 1 when $s_{i}$ and $s_{j}$ are in the same cluster and $s_{j}$ is the
parent node of $s_{i}$ in the data forwarding tree, otherwise $x_{ij}$ is set
to 0. Next, the binary variable $r_{i}$ takes value 1 if $s_{i}$ and $s_{j}$
are two CHs visited by the same UAV and $s_{j}$ is visited right after $s_{i}$
on the trajectory, or $r_{ij}=0$ if one of these above conditions is violated.

Since we aim to maximize the network lifetime, the objective function in the
MILP model must be relevant to the residual energy of SNs. Thus, we define
$n+1$ real variables $e_{1},e_{2},...,e_{n},e_{min}$ where the first $n$
variables equal to $n$ values of SN’s residual energy, respectively and the
$(n+1)^{th}$ variables that will take the minimum value among these former.

With these above variables, the problem objective and basic constraints are
given from (1) to (5).

$\begin{array}[]{l@{}lr}\textbf{Objective:$\;\;\;\;\;\;$}&\displaystyle\max
f=\alpha e_{min}+(1-\alpha)\frac{1}{N}\sum_{i\in V^{*}}e_{i}&\text{(O)}\\\
\textbf{Constraints:}&&\\\ \end{array}$

$\begin{array}[]{clr}\displaystyle\sum_{j\in
V\setminus\\{i\\}}x_{ij}+\sum_{j\in V\setminus\\{i\\}}r_{ij}=1&\forall i\in
V^{*}&\text{(B.1)}\\\ \displaystyle\sum_{j\in
V\setminus\\{i\\}}r_{ij}=\sum_{j\in V\setminus\\{i\\}}r_{ji}&\forall i\in
V^{*}&\text{(B.2)}\\\ \displaystyle\sum_{j\in V^{*}}r_{0j}\leq
m&&\text{(B.3)}\\\ x_{ij}=0&\forall i,j\in V&\\\ &\text{and
}d_{ij}>R&\text{(B.4)}\\\ x_{ij},r_{ij}\in\\{0,1\\}&\forall i,j\in
V&\text{(B.5)}\\\ e_{min}\geq 0,e_{i}\geq e_{min}&\forall i\in
V&\text{(B.6)}\\\ \end{array}$

As shown in (O), the objective of the MILP model is to maximize a weighted sum
of $e_{min}$ and $\frac{1}{N}\sum_{i\in V^{*}}e_{i}$, that are respectively
the minimum and the average residual energy of SNs. Intuitively, since the
network lifetime is equal to the time until the value of $e_{min}$ reaches 0,
we want to decrease the reduction rate of $e_{min}$ in the long-term thus our
approach is to greedily maximize the value of $e_{min}$ after each round.
However, we also need to consider the average residual energy in the network
to avoid the case of inefficient consumption at other SN with higher energy.

For example, Fig. 1 illustrates two different cluster configuration for a
network scenario where $s_{1}$ and $s_{2}$ are CHs. In this example, we assume
that all SNs start with the same energy and generate the same 1 byte of data.
Both cluster configurations yield the same value of $e_{min}$ which equal the
residual energy of the CH nodes since each of them consume the same amount of
energy to receive 1 byte and transmit 2 bytes to the UAV. However, the second
configuration is more efficient because the other nodes $s_{3}$ and $s_{4}$
consumed less energy due to shorter transmission range. Therefore, by
including the average SN energy in the objective as in (O), the better
configuration is chosen.

Figure 1: An example of two cluster configurations.

In section IV, we will do the experiment to find the optimum weight value
$\alpha$ so that the network lifetime is maximized.

Next, consider the basic constraints, constraint (B.1) means that a SN must be
either forward data to its parent node if it is not a CH
($\displaystyle\sum_{j\in V\setminus\\{i\\}}x_{ij}=1$) or else forward data to
a corresponding UAV ($\displaystyle\sum_{j\in V\setminus\\{i\\}}r_{ij}=1$).
Constraint (B.2) require that number of UAVs that arrive and leave a SN must
be the same, namely equal to 1 if it is a CH and equal to 0 otherwise. The
maximum number of UAVs is met due to (B.3). Finally, constraint (B.4)
indicates that a SN $s_{i}$ could send data to $s_{j}$ only if the Euclide
distance $d_{ij}$ between them is not exceed the maximum radio range $R$.

In the next two subsection, we present two completed MILP formulations for the
multi-UAV data collection problem by adding other constraints on SN energy and
trajectory length of UAV.

###  III-B 2-index formulation

In this formulation, we define the integer variables $B^{n}_{i}$ and
$B^{e}_{ij}$ to calculate the energy expense on each SN. In particular,
$B^{n}_{i}$ count the amount of data in bits that $s_{i}$ received from its
descendants, while $B^{e}_{ij}$ count the amount of data sent from $s_{i}$ to
$s_{j}$. The constraints related to $B^{n}_{i}$ and $B^{e}_{ij}$ are given in
(F1.1) and (F1.2).

$\begin{array}[]{clr}\displaystyle\sum_{j\in V}B^{e}_{ji}=B^{n}_{i}&\forall
i\in V^{*}&\text{(F1.1)}\\\ B^{n}_{i}+b_{i}\leq B^{e}_{ij}+M(1-x_{ij})&\forall
i\in V^{*}&\\\ &,j\in V&\text{(F1.2)}\\\ \end{array}$

Constraint (F1.1) means that the amount of data that node $s_{i}$ received
must be equal to the total amount of data sent through all of its incoming
link. In (F1.2), the big-M notation is used to formulate a disjunctive
constraint. Specifically, if $s_{j}$ is the parent node of $s_{i}$
($x_{ij}=1$) then the amount of data sent through the link $s_{i}-s_{j}$ must
not be less than the total amount that $s_{i}$ received $B^{n}_{i}$ plus the
amount generated by itself $b_{i}$. Otherwise, if $x_{ij}=0$ then the
inequality (8) becomes $B^{n}_{i}+b_{i}\leq B^{e}_{ij}+M$ and always gets true
since $M$ is chosen as a very large value.

By applying the results in (3) with the variables $B^{n}_{i}$ and
$B^{e}_{ij}$, the energy consumption at SNs can be formulated as follow:

$\begin{array}[]{cl}\displaystyle E^{r}_{i}-E_{Rx}B^{n}_{i}-\sum_{j\in
V^{*}\setminus\\{i\\}}E_{Tx}^{ij}B^{e}_{ij}&\geq e_{i}\\\ &\forall i\in
V^{*}\;\;\;\;\text{(F1.3)}\\\ \displaystyle
E^{r}_{i}-E_{Rx}B^{n}_{i}-E_{Tx}^{UAV}(B^{n}_{i}+b_{i})&\\\
\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\displaystyle\geq e_{i}-M\left(1-\sum_{j\in
V\setminus\\{i\\}}r_{ij}\right)&\forall i\in V^{*}\;\;\;\;\text{(F1.4)}\\\
\end{array}$

Here in (F1.3), the residual energy of a SN $s_{i}$ must be less or equal than
the initial value subtracted to the receiving energy ($E_{Rx}B^{n}_{i}$) and
the total transmitting energy on the outgoing links ($\sum_{j\in
V^{*}\setminus\\{i\\}}E_{Tx}^{ij}B^{e}_{ij}$). The constraint (F1.3) can be
applied for all SNs, however it only takes into account the energy consumption
for the radio communication between SNs. Then for the CHs, we also need to
consider the energy spent to transmit data to the UAV. This is done by
constraint (F1.4), which is only applied with the CH node where $\sum_{j\in
V\setminus\\{i\\}}r_{ij}=1$ by using big-M formulation.

Consider the constraint on UAV trajectory, we introduce the real variables
$L_{i}$ corresponding to each SN $s_{i}$. If $s_{i}$ is CH then $L_{i}$ is
equal to the length that its corresponding UAV (i.e. the UAV assigned with
$s_{i}$) has flown from departure until it arrive at the hovering position
above $s_{i}$. In the opposite case, the value of $L_{i}$ is set to 0.
Therefore, if $s_{i}$ and $s_{j}$ are two consecutive nodes in a UAV
trajectory (i.e. $r_{ij}=1$), the value of $L_{j}$ is equal to $L_{i}$ plus
the Euclide distance $d_{ij}$ between these two nodes. As a consequence, the
total length of an UAV trajectory is calculated as $L_{i}+d_{i0}$ where
$s_{i}$ is the last node in the trajectory (this condition is equivalent to
$r_{i0}=1$). Finally, the linear constraints related to the UAV trajectory are
given in (F1.5) and (F1.6).

$\begin{array}[]{clr}L_{i}+d_{ij}\leq L_{j}+M(1-r_{ij})&\forall i,j\in
V&\text{(F1.5)}\\\ \displaystyle L_{i}+d_{i0}r_{i0}\leq L_{max}&\forall i\in
V&\text{(F1.6)}\\\ \end{array}$

In summary, the 2-index formulation for the data collection problem with
multiple UAVs is:

$\begin{array}[]{l@{}lr}\textbf{Objective$\;\;\;\;\;\;$}&\displaystyle\max
f&\text{(F1)}\\\ \textbf{subject to }&\text{(B.1) - (B.6) and (F1.1) -
(F1.6)}&\\\ \end{array}$

###  III-C 3-index formulation

Although the above formulation is compact which only has $O(n^{2})$ variables
and constraints, however it is based on the big-M notation. It is known that
the big-M contraint typically yields a very poor lower bound in the search
trees and thus leading to a slower solving time . Thus, in this subsection, we
present a novel MILP formulation that can exclude the big-M notations.

This formulation is inspired by the multi-commodity flow formulation for the
Ring-Tree problem in [18]. Here we assume that there are $n$ flows of
commodity, each of them originates from a SN and terminates at the sink node.
We call the flow that start from a source node $s_{k}$ as the $k^{th}$
commodity flow. Then the continuous flow variable $f_{ij}^{k}\in[0,1]$ is
introduced to model the amount of the $k^{th}$ flow that is transferred from
$s_{i}$ to $s_{j}$. All flow variables must be satisfied the following
constraints:

$\begin{array}[]{clr}f_{ij}^{k}\leq x_{ij}+r_{ij}&\forall i,j,k\in
V&\text{(F2.1)}\\\ \displaystyle\sum_{j\in V^{*}}f_{0j}^{k}-\sum_{j\in
V^{*}}f_{j0}^{k}=-1&\forall k\in V^{*}&\text{(F2.2)}\\\
\displaystyle\sum_{j\in V\setminus\\{k\\}}f_{kj}^{k}-\sum_{j\in
V\setminus\\{k\\}}f_{jk}^{k}=1&\forall k\in V^{*}&\text{(F2.3)}\\\
\displaystyle\sum_{j\in V\setminus\\{i\\}}f_{ij}^{k}-\sum_{j\in
V\setminus\\{i\\}}f_{ji}^{k}=0&\forall i,k\in V^{*}&\text{(F2.4)}\\\
f_{ij}^{k}\in[0,1]&\forall i,j,k\in V&\text{(F2.5)}\\\ \end{array}$

Constraint (F2.1) means that a commodity flow can only be transferred from
$s_{i}$ to $s_{j}$ if the edge $(s_{i},s_{j})$ is either in a cluster tree
($x_{ij}=1$) or in a UAV trajectory ($r_{ij}=1$). Constraints (F2.2), (F2.3)
and (F2.4) are the flow conservation constraints at the sink node, the source
node and the transit nodes, respectively.

The flow-based model helps to eliminate the case where a solution contains
subtour in a cluster or UAV trajectory since every flows must be end at the
sink node. Moreover, these flow variables could be used to formulate other
specific constraints of our problem, namely we can calculate the energy
consumption of a SN based on the amount of flow that go in and out at that
node and calculate the UAV trajectory lengths based on the paths of some
certain flows.

First, consider the energy constraints, we need to calculate the amount of
data that a SN received and sent. This is directly related to measuring the
amount of flow that go in and out at a SN. However, we should notice that
there are two types of flow at a SN $s_{i}$: the intra-cluster flow that
originated from a SN $s_{j}$ in the same cluster with $s_{i}$ and the inter-
cluster flow that comes from a different cluster. In the real network context,
an intra-cluster flow indicates that there is a data transmission between SNs,
while an inter-cluster flow can be seen as an UAV movement between two CHs.

Figure 2: Example of intra-cluster and inter-cluster flow.

Accordingly, we define the second flow variables $g_{ij}^{k}$ to measure the
intra-cluster flows transferred between nodes. It is easy to see that a flow
transferred from $s_{i}$ to $s_{j}$ is an inter-cluster flow if and only if
$s_{i}$ and $s_{j}$ are two consecutive nodes on a UAV trajectory or when
$r_{ij}=1$ in other word (as illustrated in Fig. 2). Thus, the mathematical
form of $g_{ij}^{k}$ can be written as:

$g_{ij}^{k}=\begin{cases}f_{ij}^{k}&\text{if }r_{ij}=0\\\ 0&\text{if
}r_{ij}=1\end{cases}$ (5)

The above formula can be translated to the linear constraints as follow:

$\begin{array}[]{clr}g_{ij}^{k}\leq f_{ij}^{k}&\forall i,j,k\in
V&\text{(F2.6)}\\\ g_{ij}^{k}\geq f_{ij}^{k}-r_{ij}&\forall i,j,k\in
V&\text{(F2.7)}\\\ g_{ij}^{k}\leq 1-r_{ij}&\forall i,j,k\in V&\text{(F2.8)}\\\
g_{ij}^{k}\geq 0&\forall i,j,k\in V&\text{(F2.9)}\\\ \end{array}$

More details, when $r_{ij}=0$ then the inequality (F2.7) become
$g_{ij}^{k}\geq f_{ij}^{k}$ and by combining with (F2.6) we have
$g_{ij}^{k}=f_{ij}^{k}$. Similarly, when $r_{ij}=0$ then we have
$g_{ij}^{k}\leq 0$ from (F2.8) and $g_{ij}^{k}\geq 0$ from (F2.9), thus
$g_{ij}^{k}=0$. In both case, the value of $g_{ij}^{k}$ always satisfies (5).

With the introduction of intra-cluster flow variables $g_{ij}^{k}$, the
constraints on energy consumption of SNs are:

$\begin{array}[]{cl}\displaystyle E^{r}_{i}-E_{Rx}\sum_{j\in
V\setminus\\{i\\}}\sum_{k\in V^{*}}b_{k}g_{ji}^{k}&\\\
\displaystyle\;\;\;\;\;\;\;\;-\sum_{j\in
V\setminus\\{i\\}}E_{Tx}^{ij}\left(\sum_{k\in V^{*}}b_{k}g_{ij}^{k}\right)\geq
e_{i}&\forall i\in V^{*}\;\;\;\text{(F2.10)}\\\ \end{array}$

$\begin{array}[]{cl}\displaystyle E^{r}_{i}-(E_{Rx}+E_{Tx}^{UAV})\sum_{j\in
V\setminus\\{i\\}}\sum_{k\in
V^{*}\setminus\\{i\\}}&b_{k}(g_{ji}^{k}-g_{ij}^{k})\\\ \displaystyle-
E_{Tx}^{UAV}b_{i}\sum_{j\in V\setminus\\{i\\}}r_{ij}\geq e_{i}&\forall i\in
V^{*}\;\;\;\text{(F2.11)}\\\ \end{array}$

Next, the constraints on UAV trajectory lengths are given in (F2.12)

$\begin{array}[]{c}\displaystyle d_{0k}\sum_{j\in
V\setminus\\{i\\}}r_{ij}+d_{ij}\sum_{j\in
V\setminus\\{i\\}}(f_{ij}^{k}-g_{ij}^{k})\leq L_{max}\end{array}$
$\begin{array}[]{r}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\forall
k\in V^{*}\;\;\;\text{(F2.12)}\end{array}$

To sum up, the proposed 3-index formulation is as follow:

$\begin{array}[]{l@{}lr}\textbf{Objective$\;\;\;\;\;\;$}&\displaystyle\max
f&\text{(F2)}\\\ \textbf{subject to }&\text{(B.1) - (B.6) and (F2.1) -
(F2.12)}&\\\ \end{array}$

###  III-D Formulation analysis

Both presented formulation employ $2n^{2}+n+1$ base variables including
$n^{2}$ variables $x_{ij}$, $n^{2}$ variables $r_{ij}$ and $n+1$ variables
$e_{min},e_{1},...,e_{n}$. In the 2-index formulation, there are an additional
of $n^{2}+2n$ variables including $Bn_{i},Be_{ij}$ and $L_{i}$. Thus there are
in total $3n^{2}+3n+1$ variables in this formulation. For the 3-index
formulation, we introduce $2n^{3}$ flow variables $f^{k}_{ij}$ and
$g^{k}_{ij}$. However, since we have $x_{ij}=0$ and $g^{k}_{ij}=0$ for all
$i,j$ such that $d_{ij}>R$ (i.e. $s_{j}$ is out of transmission range of
$s_{i}$) then we can reformulate the 3-index formulation so that there are
only $n\delta$ variables $x_{ij}$ and $n^{2}\delta$ variables $g^{k}_{ij}$
where $\delta$ is the maximum number of neighbor nodes of a SN. In general,
the 3-index formulation has higher number of $O(n^{3})$ variables and
constraints compared to $O(n^{2})$ in the 2-index formulation, however we will
show in the experiment that the solving time of the 3-index formulation is
significantly faster than the other one because it can exclude the big-$M$
notation.

###  III-E Heuristic model

Through the experiment, we find that our best MILP model can only deal with
the problem of moderate size, namely up to 50 SNs thus it is insufficient to
apply in a large-scale network of several hundred or thoudsands of nodes. Thus
we propose a simple heuristic to reduce the number of variables in the MILP
model. Since our objective is to prolong the network lifetime, the selected
CHs in a round should be the nodes with higher residual energy. Therefore, we
restrict the set of potential nodes to become CHs to the set of $P$ percent of
nodes with the highest energy in the network $N_{P}$. As a result, the number
of variables $r_{ij}$ can be reduced from $n^{2}$ to $|N_{P}|^{2}$ because
$r_{ij}=0$ if either i or j is not in the set of potential CHs $N_{P}$ and the
number of variables $f^{k}_{ij}$ can also be reduced from $n^{3}$ to
$|N_{P}|^{3}$.

The first competitor is the Convex-Hull based protocol (CHP) in [10], which
also handled the case of multiple UAVs. This protocol firstly divide the
network into multiple sectors with the same percentage of energy decrease and
each mobile collector (MC) is associated with an unique sector. An UAV
trajectory is constructed inside its sector at initial by a geometrical method
and then a simulated-annealing adjustment is adopted to improve the
trajectories in term of total network energy consumption.

##  IV Experiments

In this section, we first compare internally between two proposed MILP models
and the heuristic implementation regarding the computation time and
optimality. Later we evaluate the network lifetime performance when applying
the heuristic model in comparison with an existing work. In the experiments,
the network size is varied from 15 to 100 nodes. The communication range of a
sensor node is 40m, and the flying altitude $H$ of each UAV is 20m. By
default, each SN $s_{i}$ generates a data packet of a fixed size of 100Kb per
round. Each SN’s initial energy is randomly selected from the interval [20J,
40J], thus reflecting heterogeneous WSNs. Table I sums up the experiment
settings and parameters.

All MILP models were solved by using the commercial MIP solver CPLEX in
version 20.10. The experiments were run on a computer with an Intel(R)
Core(TM) i7-4770 CPU @ 3.40GHz and 8 GB of RAM. For each test instance, the
relative MIP gap tolerance is set to 0.1%. We also set a solving time limit of
30 minutes in CPU time. Then if the running time for a test instance exceeds
the predefined time limit, we terminate the solving process and record the
best gap achieved. Other configuration parameters of the CPLEX solver are left
as default.

TABLE I: Default experiment settings Network | Area | 15/20/30 SNs: 50m $\times$ 50m
---|---|---
50 SNs: 75m $\times$ 75m
| 100 SNs: 100m $\times$ 100m
Sink location | Center of the square area
Number of SNs | 15 - 100
| Number of UAVs | 2
| Max. trajectory length | 100m
| Packet size | 100Kb
| SN energy | Random in [20J, 40J]
Radio | Communication range | 40m
$E_{elec}$ | 50 nJ/bit
$\epsilon_{fs}$ | 10 pJ/bit/m2
| $\epsilon_{mp}$ | 0.0013 pJ/bit/m4

###  IV-A Impact of the weight parameter $\alpha$

Firstly, we examine the most efficient value of the weight factor $\alpha$ in
the objective function that would induce the longest network lifetime. In
particular, we experiment on a random WSN with 30 SNs and the value $\alpha$
is varied from 0.1 to 1. The network keep operating until the first SN is
exhausted and we apply the 2-index model to obtain the network configuration
in each round. Besides, we use a smaller data packet size of 10Kb for all SNs
to differentiate between results more clearly. Fig 3 depicts the experiment
results of network lifetime with different value of $\alpha$. As can be seen,
the value $\alpha$ equal to 0.6 achieved the longest lifetime, namely
approximate 2000 rounds. Accordingly, we will used the value $\alpha=0.6$ in
the next experiments.

Figure 3: Impact of the weight $\alpha$ to the network lifetime

###  IV-B Comparison between models

In this second experiment, we will compare the performance of two proposed
MILP models and the heuristic implementation in terms of solving time and
optimal gap. Here, we test each model on several network instances where the
number of SNs and there are 20 network instances randomly generated for each
network size. In this experiment, we also consider two different energy
settings. In particular, for each network size, there are only 10 instances
that take the default energy setting as presented above. With the 10 remaining
network instances, each SN will take a random energy value from a more narrow
interval [10J, 12J]. This second energy setting reflects the situation in a
later round when the SN energies are gradually balanced through multiple
previous rounds.

We provide the comprehensive results in Table II. In this table, the solving
time and the best relative gap to the optimal solution, including the average
and the worst case value, are respectively presented in the column ”time” and
”gap”. Besides, the column ”$P_{opt}$” shows the percentage of instances that
can be solved optimally over 20 instances. Especially, a cell in column ”gap”
or ”$P_{opt}$” is left blank in case that all test instances are solved
optimally.

According to the results,

TABLE II: Comparison between models N | Model | $P_{opt}$ | Time (s) | Gap (%)
---|---|---|---|---
avg | max | avg | max
15 | 2-index |  | 0.50 | 2.05 |  |
3-index |  | 2.83 | 7.75 |  |
Heuristic 50% |  | 0.13 | 0.20 |  |
Heuristic 20% |  | 0.02 | 0.14 |  |
30 | 2-index |  | 5.58 | 20.23 |  |
3-index | 45% | 1026.02 | 1800.00 | 2.07 | 9.24
Heuristic 50% |  | 2.63 | 7.17 |  |
Heuristic 20% |  | 0.37 | 0.67 |  |
50 | 2-index | 95% | 256.66 | 1800.00 | 0.05 | 0.10
3-index | 0% | 1800.00 | 1800.00 |  |
Heuristic 50% |  | 38.62 | 122.77 |  |
Heuristic 20% | 95% | 5.82 | 12.17 | 0.04 | 0.11
100 | 2-index |  |  |  |  |
3-index |  |  |  |  |
Heuristic 50% |  |  |  |  |
Heuristic 20% |  | 126.61 | 360.09 |  |

###  IV-C Performance of the heuristic model

To illustrate the effectiveness of our heuristic MILP approach, we compare it
with the Convex-Hull-based protocol (CHP) in [10]. This protocol firstly
divides the network into multiple sectors with the same percentage of energy
decrease, and each mobile collector (UAV) is associated with a unique sector.
A UAV trajectory is constructed inside its sector at initial by a geometrical
method. Then a simulated-annealing adjustment is adopted to improve the
trajectories in terms of total network energy consumption.

Fig. 4 presents the minimum residual energy after each round when applying our
heuristic model and the CHP in a network with 100 SNs. As can be seen, our
algorithm has a lower rate of energy reduction. As a result, the network
lifetime when applying the heuristic model is 150 rounds, higher than 36%
compare to 110 rounds as in the CHP.

Figure 4: Network lifetime in multi UAV scenario

##  V Conclusion

## Acknowledgment

This research is funded by Vingroup Joint Stock Company and supported by the
Domestic Master/PhD Scholarship Programme of Vingroup Innovation Foundation
(VINIF), Vingroup Big Data Institute (VINBIGDATA); funded by Vietnam National
Foundation for Science and Technology Development (NAFOSTED) under grant
number 102.02-2017.316.

## References

  * [1] Y. Gu, F. Ren, Y. Ji, and J. Li.  The evolution of sink mobility management in wireless sensor networks: A survey.  IEEE Communications Surveys Tutorials, 18(1):507–524, 2016.
  * [2] M. Ma, Y. Yang, and M. Zhao.  Tour planning for mobile data-gathering mechanisms in wireless sensor networks.  IEEE Transactions on Vehicular Technology, 62(4):1472–1483, 2013\.
  * [3] X. Xu, Y. Zeng, Y. L. Guan, and R. Zhang.  Overcoming endurance issue: Uav-enabled communications with proactive caching.  IEEE Journal on Selected Areas in Communications, 36(6):1231–1244, 2018.
  * [4] C. Zhan, Y. Zeng, and R. Zhang.  Energy-efficient data collection in uav enabled wireless sensor network.  IEEE Wireless Communications Letters, 7(3):328–331, 2018.
  * [5] G. Xing, T. Wang, Z. Xie, and W. Jia.  Rendezvous planning in wireless sensor networks with mobile elements.  IEEE Transactions on Mobile Computing, 7(12):1430–1443, 2008.
  * [6] Mianxiong Dong, Kaoru Ota, Man Lin, Zunyi Tang, Suguo Du, and Haojin Zhu.  Uav-assisted data gathering in wireless sensor networks.  The Journal of Supercomputing, 70(3):1142–1155, Dec 2014.
  * [7] S. Alfattani, W. Jaafar, H. Yanikomeroglu, and A. Yongacoglu.  Multi-uav data collection framework for wireless sensor networks.  In 2019 IEEE Global Communications Conference (GLOBECOM), pages 1–6, 2019.
  * [8] Qiyue Wu, Peng Sun, and Azzedine Boukerche.  Unmanned aerial vehicle-assisted energy-efficient data collection scheme for sustainable wireless sensor networks.  Computer Networks, 165:106927, 2019.
  * [9] Amar Kaswan, Kumar Nitesh, and Prasanta K. Jana.  Energy efficient path selection for mobile sink and data gathering in wireless sensor networks.  AEU - International Journal of Electronics and Communications, 73:110 – 118, 2017.
  * [10] D. Ebrahimi, S. Sharafeddine, P. Ho, and C. Assi.  Data collection in wireless sensor networks using uav and compressive data gathering.  In 2018 IEEE Global Communications Conference (GLOBECOM), pages 1–7, 2018.
  * [11] Charalampos Konstantopoulos, Nikolaos Vathis, Grammati Pantziou, and Damianos Gavalas.  Employing mobile elements for delay-constrained data gathering in wsns.  Computer Networks, 135:108 – 131, 2018.
  * [12] M. Ma and Y. Yang.  Data gathering in wireless sensor networks with mobile collectors.  In 2008 IEEE International Symposium on Parallel and Distributed Processing, pages 1–9, 2008.
  * [13] M. Zhao and Y. Yang.  Bounded relay hop mobile data gathering in wireless sensor networks.  IEEE Transactions on Computers, 61(2):265–277, 2012.
  * [14] C. Zhan, Y. Zeng, and R. Zhang.  Energy-efficient data collection in uav enabled wireless sensor network.  IEEE Wireless Communications Letters, 7(3):328–331, 2018.
  * [15] Raj Anwit, Abhinav Tomar, and Prasanta K. Jana.  Tour planning for multiple mobile sinks in wireless sensor networks: A shark smell optimization approach.  Applied Soft Computing, 97:106802, 2020.
  * [16] Rumpa Dasgupta and Seokhoon Yoon.  Energy-efficient deadline-aware data-gathering scheme using multiple mobile data collectors.  Sensors, 17(4), 2017.
  * [17] W. R. Heinzelman, A. Chandrakasan, and H. Balakrishnan.  Energy-efficient communication protocol for wireless microsensor networks.  In Proceedings of the 33rd Annual Hawaii International Conference on System Sciences, pages 10 pp. vol.2–, 2000.
  * [18] Alessandro Hill and Silvia Schwarze.  Exact algorithms for bi-objective ring tree problems with reliability measures.  Comput. Oper. Res., 94(C):38–51, June 2018.


# Inhomogeneous Floquet thermalization

Soumya Bera  Department of Physics, Indian Institute of Technology Bombay,
Mumbai 400076, India  Ishita Modak  Department of Physics, Indian Institute of
Technology Bombay, Mumbai 400076, India  Roderich Moessner  Max-Planck-
Institut für Physik komplexer Systeme, 01187 Dresden, Germany

###### Abstract

How a closed system thermalizes, especially in the absence of global
conservation laws but in the presence of disorder and interactions, is one of
the central questions in non-equilibrium statistical mechanics. We explore
this for a disordered, periodically driven Ising chain. Our numerical results
reveal inhomogeneous thermalization leading to a distribution of
thermalization timescales within a single disordered sample, which we encode
via a distribution of effective local temperatures. Using this, we find an
excellent collapse without any fitting parameters of the local relaxation
dynamics for the entire range of disorder values in the ergodic regime when
adapting the disorder-averaged diagonal entanglement entropy as internal
‘time’ of the system. This approach evidences a remarkably uniform
parametrization of the dynamical many-body evolution of local temperature
within the otherwise highly heterogeneous ergodic regime, independent of the
strength of the disorder.

#### Introduction:

Spatially heterogeneous relaxation dynamics towards equilibrium is a hallmark
of nonergodicity, being found in paradigmatic settings of glasses and jammed
systems [1]. Such dynamical heterogeneity, e.g., evidenced in the coexistence
of different relaxation timescales, can arise from spatial variations
associated with the presence of metastable states. This can lead to global
nonexponential decay of correlation functions in time despite local
exponential decay rates. Alternatively, relaxation processes can be inherently
complex also, contributing to local nonexponential decay [2]. In ergodic
systems, such as supercooled liquids, even regions of slower relaxation
eventually thermalize [3].

In quantum systems undergoing unitary dynamics, the nature of thermalization,
or its absence, has been a focal point of research in recent years [4, 5, 6].
Thermalization in ergodic systems occurs, loosely speaking, as sub-system
density matrices evolve to a thermal state, with the remainder of the system
effectively acting as its bath. The presence of disorder may impede
thermalization, for instance, by the emergence of a set of quasi-local
integrals of motion, resulting in emergent effective integrability and non-
ergodicity, a phenomenon dubbed many-body localization (MBL) [7, 8, 4, 9, 10,
11, 12]. In this parameter regime, spatial and dynamical heterogeneity has
been observed in local entanglement measures [13, 14] as demonstrated in
models such as the XXZ spin chain. Indeed, entanglement can exhibit
substantial spatial heterogeneity in this regime, as indicated by the
subvolume scaling of the standard deviation of the cut-to-cut entanglement
entropy [15]. Because of the difficulty of studying the real-time dynamics of
disordered interacting quantum systems, a definite consensus on the nature of
this regime has been slow to emerge, e.g., Refs. [16, 17, 18, 19, 20, 21, 22,
23, 24, 25, 26, 27, 28, 29].

Here, we show that even the weakly disordered ergodic regime can exhibit
considerable spatial structure, which we investigate in detail. We focus on a
nonintegrable driven disordered Ising chain without global conservation laws.
While even weak disorder tends to slow down relaxation, eventual
thermalization can be remarkably robust [16]. We analyze thermalization after
a sudden quench via the time evolution of the local subsystem (inverse)
temperature, $\beta_{j}(n)$, where $j$ denotes a bond involving two sites
after $n$ the time steps. The $\beta_{j}(n)$ evolution reveals apparent
glassiness in the sense of a distribution of local relaxation timescales
within a single disorder configuration. Indeed, the time evolution of the
spatial and disorder averaged $\beta_{j}(n)$ shows well-developed
nonexponential behavior in finite-size numerical simulation. With increasing
disorder, we further observe mixed dynamics, i.e., locally nonexponential
decay of the $\beta_{j}(n)$ accompanied by regions in space with exponential
relaxation in time towards ‘infinite temperature’.

The broad distribution of thermalization times, even at a single disorder
value, suggests that relying on a single time scale may not be sufficient.
Relatedly, it is challenging to compare the time evolution between different
disorder realizations and, further, strengths. Here, we propose a unified
description employing the ensemble-averaged diagonal entanglement entropy
${{\mathscr{S}}}_{\text{d}}(t)$ to monitor the dynamical evolution of
$\beta_{j}(n)$ across the entire ergodic regime. This extends the framework
introduced by Evers  _et al._ [28] for non-driven systems, which used the
average entanglement entropy to track particle density imbalance decay. We
find a straightforward data collapse of the evolution of average inverse
temperature for several disorder values, offering a homogeneous perspective on
thermalization dynamics despite its inherent heterogeneity. While previous
work [28] required extra fitting parameters to obtain a scaling collapse, our
intrinsic parametrization of time obviates this need.

#### Driven Ising Model:

The time-dependent Hamiltonian of this periodically driven system is defined
as,

$\displaystyle H(t)=\begin{cases}2H_{x}&\text{if }0<t<\frac{T}{2}\\\
2H_{z}&\text{if }\frac{T}{2}<t<T,\end{cases}$ $\displaystyle
H_{x}=\sum_{i=1}^{L}g\Gamma\sigma_{i}^{x},$ $\displaystyle
H_{z}=J\sum_{i=1}^{L-1}\sigma_{i}^{z}\sigma_{i+1}^{z}+\sum_{i=1}^{L}(h+g\sqrt{1-\Gamma^{2}}G_{i})\sigma_{i}^{z},$
(1)

$\sigma_{i}^{x}$ and $\sigma_{i}^{z}$ are the Pauli matrices on site $i$. We
follow the standard parametrization of the model, see Refs. [30, 31, 32]. The
interaction strength $J=1$ and $g,h$ and $T$ are 0.9045, 0.8090 and 0.8,
respectively. Such a parametrization is motivated by the clean static model,
where strong thermalization is obtained with these parameter values for system
sizes readily accessible in exact diagonalization studies [33]. The
longitudinal field is disordered and chosen from a Gaussian distributed random
variable $G_{i}$ with zero mean and unit variance. The $\Gamma$ controls the
disorder strength, and the model is believed to have an MBL transition at
$\Gamma\simeq 0.3$ [30, 31]. The stroboscopic time evolution is performed with
Floquet operator $U_{\text{F}}(T)=e^{-iH_{x}T/2}e^{-iH_{z}T/2}$ using a
Hadamard transformation, see Ref. [31] for further details. The initial state
for these calculations is the Néel state $|1010\ldots\rangle$.

Figure 1: Stretched exponential decay of the average local temperature
$\overline{\beta}(n)$ with Floquet time step $n$. The line indicates a
stretched exponent fit: $A\exp[-(t/\tau)^{\alpha})]$. The corresponding values
of the thermalization time $\tau(L)$ and the stretched exponent $\alpha(L)$
are provided in the legend for $L=20$. Lower panel: Spatiotemporal
inhomogeneity of the evolution of local bond temperature $\beta_{j}(n)$ as
defined in Eq. (2) for different values of $\Gamma={0.7,0.6,0.5}$ as a
function of $\log(n)$ for a typical disorder configuration. The white coloring
represents infinite temperature ($\beta_{j}{=}0$).  Figure 2: Time evolution
of local inverse temperature $|\beta_{j}(n)|$ at different locations, $j$, for
several disorder values $\Gamma=0.6,0.5,0.4,0.45$ for $L=20$. Log-linear scale
highlights exponential decay at large $\Gamma$ for a typical disorder
configuration. With decreasing $\Gamma$, the thermalization time increases.
For the small value of $\Gamma=0.45$, the curves appear to fall into two
groups characterized by rapid and slow decay rates.

#### Local temperature:

The local temperature $\beta_{j}(n)$ is defined for each bond by minimizing
the Frobenius norm distance between the canonical density matrix
$e^{-\beta_{j}\mathscr{H}_{\text{b}}}$, and the sub-system density matrix
$\rho_{j}^{\text{b}}$,

$\displaystyle\text{min}_{\beta_{j}}\,\left\lvert
e^{-\beta_{j}\mathscr{H}_{\text{b}}}-\rho_{j}^{\text{b}}(n)\right\rvert\ ,$
(2)

with
$\mathscr{H}_{\text{b}}(j)=\sigma^{z}_{j}\sigma^{z}_{j+1}+h_{j}/2\,\sigma^{z}_{j}\otimes\mathds{1}_{2\times
2}+h_{j+1}/2\,\mathds{1}_{2\times 2}\otimes\sigma^{z}_{j}+\\{z\leftrightarrow
x,h_{j}\leftrightarrow g_{j}\\}$ is the local bond Hamiltonian with
$h_{j},g_{j}$ being the local fields, and
$\rho^{\text{b}}_{j}(n)=\text{Tr}_{L-\\{j,j+1\\}}\,|\psi(n)\rangle\langle\psi(n)|$
is the reduced density matrix for that bond. The norm is defined as $\lvert
A\rvert\equiv\text{tr}(\sqrt{A^{\dagger}A})$. We expand the definition of
$\beta$ from Ref. [34] to the time domain. This approach includes the time-
evolved wavefunction’s structure through the sub-system density matrix. It is
important to note that the precise value of $\beta_{j}(n)$ is contingent upon
the chosen definition of norm, as examined in detail in Ref. [34].

#### (Non-)Exponential Heating:

An interacting driven system heats to an infinite temperature featureless
state [35, 36]. The heating rate $\tau$ generally depends exponentially on the
drive frequency, $\tau\propto\exp(\omega/J)$ for $\omega/J\gg 1$, where $J$ is
a microscopic energy scale [37, 38, 39, 40]. In the presence of disorder, it
has been seen that the (average) heating slows down considerably [41, 42, 43,
44, 31]. As shown in Fig. 1, the time evolution of the disorder and spatially
averaged (indicated by the overline) $\overline{\beta}(n)$ indeed exhibits a
stretched exponential decay, with the decay exponent inversely correlated to
the strength of the disorder, similar to the decay of the correlation function
as reported earlier in Refs. [31, 27, 45]. A pronounced finite size effect is
also observed and shown in the Appendix.

The lower panel of Fig. 1 highlights the spatially inhomogeneous evolution of
$\beta_{j}(n)$ for a given disorder configuration for several $\Gamma$ values.
The upper panel of Fig. 2 shows $\beta_{j}(n)$ at different locations for a
single disorder configuration. For weak disorder, $\Gamma=0.7$, thermalisation
is exponential everywhere, $e^{-n/\tau_{j}}$, however with spatially varying
$\tau_{j}$. The $|\beta_{j}(n)|$ decays to ${\sim}10^{-3}$ for $L=20$, and is
expected to vanish in the thermodynamic limit. With increasing disorder, the
variation in $\tau_{j}$ increases, and for even stronger disorder
$\Gamma=0.45$, the heating becomes slow, possibly as a stretched exponential
(see Fig. 2(c)). Even there, a few subsystems still exhibit fast
thermalization, i.e., an exponential decay of $|\beta_{j}(n)|$ in time, $n$.
This distribution of thermalization time scale is reflected in a stretched
exponential decay, $e^{-(t/\tau)^{\alpha}}=\int duP(u)e^{(-t/u)}$, of the
average $\overline{\beta}(n)$ as highlighted in Fig. 1 with a timescale that
increases, and an exponent $\alpha$ which decreases, with increasing disorder.

#### Thermalization time:

Fitting each $\beta_{j}(n)$ trace is impractical because of fluctuations in
the data and uncertainties associated with the fit. We instead extract a local
decay time $\tau_{j}$ via

$\displaystyle\tau_{j}\coloneqq\frac{\int_{0}^{T}n|\beta_{j}(n)|dn}{\int_{0}^{T}|\beta_{j}(n)|dn}.$
(3)

Concretely, the pure exponential $|\beta_{j}(n)|\propto
e^{-n/\tau_{j}^{\prime}}$ yields
$\tau_{j}=\tau_{j}^{\prime}+T/(1+e^{T/\tau_{j}^{\prime}})\rightarrow\tau_{j}^{\prime}$
for sufficiently large number of Floquet time steps $T$. Indeed, for the
larger $\Gamma=0.7$, the individual samples show exponential decay within the
simulation time window. With decreasing $\Gamma<0.5$, not all the traces of
$\beta_{j}(n)$ show pure exponential decay; instead, there is mixing of both
exponential and stretched exponential decays (Fig. 2, right panel). This
affects the calculation of the decay time as defined in Eq. (3), and we
refrain from doing this analysis for smaller values of the $\Gamma$, i.e.,
larger disorder values close to the putative MBL transition. In this regime,
the simple way of describing the exponential heating dynamics using the Fermi-
Golden rule [40, 46, 47] is probably inapplicable, and one might need to go
beyond this perturbative treatment.

Figure 3: Probability distribution of the local thermalization time $\tau_{j}$
for system sizes $L=18,20,22$, and disorder values $\Gamma=0.7,0.6,0.5$. The
distribution is broad with an exponential tail. For smaller $\Gamma=0.5$, the
distribution shows strong finite-size effects. Typically, $\sim
2.5{\times}10^{3}$ disorder configurations are used for these distributions.

Figure 3 shows the $\tau_{j}$ distribution in the ergodic phase for
$\Gamma=0.7,0.6,0.5$. A pronounced exponential tail is observed for all the
disorder values, with a plateau forming at smaller $\Gamma$, revealing a broad
distribution of timescales. At $\Gamma=0.5$, significant finite size and time
effects are evident (see also Appendix), suggesting the need for larger system
sizes $L\gtrsim 22$ to observe thermalization in all parts of the sample.

For the largest system $L=22$, the $\mathscr{P}(\tau_{j})$ shifts towards the
origin compared to smaller system sizes, indicating a slow flow of the full
distribution towards shorter thermalization time.

#### Density matrix evolution:

Having established a relatively broad distribution of $\tau_{j}$, resulting in
a stretched exponential decay of the spatially averaged $\beta_{j}(n)$ (Fig.
1), we now provide an entanglement perspective to the thermalization process.
Figure 4 shows a typical evolution of the half-system density matrix
$|\rho_{nn^{\prime}}|/\text{max}(\rho_{nn^{\prime}})$ for different evolution
times, $n$, for $\Gamma=0.7$, which becomes thermal at these times. This
behavior is representative of all samples reaching thermal equilibrium.

#### Diagonal entanglement entropy:

The entanglement entropy is defined as,
${{\mathscr{S}}}_{\text{E}}=-\text{Tr}\left(\rho^{\text{A}}\log\rho^{\text{A}}\right),$
where $\rho^{\text{A}}$ is the reduced half-system density matrix. At long
time, for the ergodic system, ${{\mathscr{S}}}_{\text{E}}$ reaches the Page
value i.e., ${{\mathscr{S}}}_{\text{E}}=L/2\log(2)-1/2$ [48]. In this regime,
the diagonal elements $\rho_{kk}^{A}$ dominate; these scale as
$\rho_{kk}^{A}\propto 1/\sqrt{\mathscr{D}}$, while the off-diagonal terms are
suppressed as $\sqrt{\mathscr{D}^{1/2}}$ [49], where $\mathscr{D}=2^{L}$ is
the Fock space dimension. The diagonal entropy in this basis is

$\displaystyle\mathscr{S}_{\text{d}}=-\sum_{k}\rho_{kk}^{A}\log\rho_{kk}^{A}.$
(4)

For a Haar random state, $\rho_{kk}^{A}\propto 2^{-L/2}$, therefore
$\mathscr{S}_{\text{d}}\propto(L/2)\log(2)$, the volume law scaling for
ergodic systems. Further, the diagonal entropy can be expressed in terms of
the participation entropy $\mathscr{S}_{\mathrm{P}}=2\mathscr{S}_{\text{d}}$
(see Appendix for details), where
${{\mathscr{S}}}_{\text{P}}(n)=-\sum_{j=1}^{\mathscr{D}}p_{j}(n)\log
p_{j}(n)$, with $p_{j}(n)=|\langle j|\psi(n)\rangle|^{2}$ is the probability
of occupation of each spin basis state $|j\rangle$, and $|\psi(n)\rangle$ is
the time evolved wavefunction.

Therefore, ${{\mathscr{S}}}_{\text{d}}$ is an alternative measure of
delocalization in configuration space. Indeed, Ref. [50] showed that for a
pure state with few random non-zero elements relative to the dimension of the
space, the scaling of ${{\mathscr{S}}}_{\text{d}}$ with subsystem volume $L/2$
is exactly given by that of its participation entropy
${{\mathscr{S}}}_{\text{P}}$.

Figure 4: Time evolution of the half-system density matrix
$\rho^{\text{A}}_{kk^{\prime}}$ (normalized such that its maximum value is
unity) for $\Gamma=0.7$, $L=12$ and for a typical disorder configuration. As
expected for a featureless random state, the diagonal elements dominate over
time as the system heats up. Additionally, the diagonal entropy reaches
maximum
$\tilde{\mathscr{S}}_{\text{d}}=\mathscr{S}_{\text{d}}/(L/2)\simeq\log(2)$, as
indicated in the figure.

#### Synchronized dynamics:

The upper panel of Fig. 5 shows the time evolution of the disorder averaged
$x(n)=-\overline{\log\beta_{\text{typ}}}(n)$ with Floquet cycle $n$ for
several values of disorder $\Gamma=0.7\ldots 0.3$. The typ. denotes the median
value of $\beta_{j}(n)$ across the bonds $j$. For intermediate disorder
values, $x(n)$ shows logarithmically slow propagation (see, e.g.,
$\Gamma=0.45$); however, with increasing time, it accelerates as seen by the
leftward bending of the curve to eventual thermalization. Such bending happens
at progressively higher $n$ with increasing disorder, and to observe this at
an even larger disorder, larger $L$ is necessary.

Most strikingly, when ${{\mathscr{S}}}_{\text{d}}$ is adapted as an ensemble
average internal time of the system, we observe an excellent collapse of the
mean $x(n)$ for several disorder values as seen in Fig. 5. This collapse,
requiring no fitting parameters, implies that the diagonal entropy
${{\mathscr{S}}}_{\text{d}}$ faithfully describes central aspects of the
thermalization of the closed system, such as the time to heat up to infinite
temperature at finite disorder values. Concretely, disorder slows down entropy
production, thus delaying thermalization. Once the simulation time is
parametrized by the entropy itself, the universal nature of the dynamics is
revealed.

Figure 5: Time evolution of the disorder averaged typical inverse local
temperature $x(n)=-\overline{\log\beta_{\text{typ}}}(n)$ for several disorder
values until the putative transition point $\Gamma_{\text{c}}\simeq 0.3$ and
for $L=20$. The initial slow growth gives way to faster thermalization at
later $n$. This feature has yet to manifest for smaller $\Gamma$. Below: The
collapse of the evolution when diagonal entropy ${{\mathscr{S}}}_{\text{d}}$
is adapted as an implicit measure of time, $n$, for the above figure. Inset:
The same data, plotted against the ${{\mathscr{S}}}_{\text{P}}$, normalized
with the Fock space dimension $\mathscr{D}$.

#### Conclusions:

For the thermalization of a disordered interacting Floquet system, we analyze
the dynamics of sub-system temperature $\beta_{j}(n)$ in the ergodic regime.
Generically, some blocks heat up faster than others, but all blocks eventually
thermalize, leading to a broad distribution of thermalization times. Blocks
with long thermalization times are not particularly _rare_ and may exhibit
either long exponential decay time constants or even nonexponential decay.
This distribution of time scales leads to a slower decay of the spatial and
disordered average inverse temperature $\beta_{j}(n)$, resulting in
nonexponential heating over time in the ergodic phase, resembling the
relaxation dynamics of classical glassy liquids. Exploring the connection
between inhomogeneous thermalization dynamics and the avalanche mechanism or
many-body resonances [51, 52, 53, 54, 55, 25, 26, 27, 56], which predict
stretched exponential decay of correlation functions, is an obvious avenue for
further study.

Identifying diagonal entropy ${{\mathscr{S}}}_{\text{d}}$ as an internal
system time allows a collapse of the thermalization dynamics across all the
disorder values in the ergodic phase, revealing a remarkable, albeit hidden,
homogeneity. The thermalization process involves shrinking off-diagonal matrix
elements of the density matrix. In this sense, the approximation of time with
diagonal entropy ${{\mathscr{S}}}_{\text{d}}$ measures the heating rate along
with the Fock space delocalization.

The prediction of any dynamical exponent is challenging due to limits imposed
by finite time and system sizes, as copiously noted in Hamiltonian models [57,
16, 17, 58, 59, 60, 28], disordered Floquet models [31, 27], and even in clean
models, where $L\gtrsim 22$ is often necessary to observe heating towards
infinite temperature [47]. Indeed, there is substantial variation in stretched
exponents with increasing system sizes, particularly evident in two-point
correlators, and their fate even in the ergodic phase in the asymptotic limit
[31, 45] is at this point unclear. We note that the data collapse we observe
is largely independent of system size and thus appears less afflicted by
finite-size effects.

#### Acknowledgements:

SB would like to thank F. Evers for several insightful discussions and
collaboration on related topics. We would also like to thank C. Artiaco and M.
Haque for several discussions. SB acknowledges support from MPG for funding
through the Max Planck Partner Group at IIT Bombay and also thanks the MPI-
PKS, Dresden computing cluster, where the calculation is performed. This work
was in part supported by the Deutsche Forschungsgemeinschaft under grant
cluster of excellence ct.qmat (EXC 2147, project-id 390858490). IM
acknowledges financial support from the Prime Minister’s Research Fellows
(PMRF) scheme offered by the Ministry of Education, Government of India.

## References

  * Berthier and Biroli [2011] L. Berthier and G. Biroli, Theoretical perspective on the glass transition and amorphous materials, Rev. Mod. Phys. 83, 587 (2011).
  * Glotzer  _et al._ [1998] S. C. Glotzer, N. Jan, T. Lookman, A. B. MacIsaac, and P. H. Poole, Dynamical heterogeneity in the ising spin glass, Phys. Rev. E 57, 7350 (1998).
  * Ediger [2000] M. D. Ediger, Spatially heterogeneous dynamics in supercooled liquids, Annu. Rev. Phys. Chem. 51, 99 (2000).
  * Nandkishore and Huse [2015] R. Nandkishore and D. A. Huse, Many-body localization and thermalization in quantum statistical mechanics, Ann. Rev. Condens. Matter Phys. 6, 15 (2015).
  * 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, Adv. Phys. 65, 239 (2016).
  * Ueda [2020] M. Ueda, Quantum equilibration, thermalization and prethermalization in ultracold atoms, Nature Reviews Physics 2, 669 (2020).
  * Gornyi  _et al._ [2005] I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Interacting electrons in disordered wires: Anderson localization and low-$t$ transport, Phys. Rev. Lett. 95, 206603 (2005).
  * Basko  _et al._ [2006] D. M. Basko, I. L. Aleiner, and B. L. Altshuler, Metal insulator transition in a weakly interacting many electron system with localized single particle states, Ann. Phys. 321, 1126 (2006).
  * Bera  _et al._ [2015] 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).
  * Abanin  _et al._ [2019] D. A. Abanin, E. Altman, I. Bloch, and M. Serbyn, Colloquium: Many-body localization, thermalization, and entanglement, Rev. Mod. Phys. 91, 021001 (2019).
  * Alet and Laflorencie [2018] F. Alet and N. Laflorencie, Many-body localization: An introduction and selected topics, C. R. Phys. (2018).
  * Sierant  _et al._ [2024] P. Sierant, M. Lewenstein, A. Scardicchio, L. Vidmar, and J. Zakrzewski, Many-body localization in the age of classical computing (2024), arXiv:2403.07111 [cond-mat.dis-nn] .
  * Bera and Lakshminarayan [2016] S. Bera and A. Lakshminarayan, Local entanglement structure across a many-body localization transition, Phys. Rev. B 93, 134204 (2016).
  * Artiaco  _et al._ [2022] C. Artiaco, F. Balducci, M. Heyl, A. Russomanno, and A. Scardicchio, Spatiotemporal heterogeneity of entanglement in many-body localized systems, Phys. Rev. B 105, 184202 (2022).
  * Khemani  _et al._ [2017] 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).
  * Weiner  _et al._ [2019] 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).
  * Panda  _et al._ [2020] 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).
  * Kiefer-Emmanouilidis  _et al._ [2020] 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).
  * Luitz and Lev [2020] D. J. Luitz and Y. B. Lev, Absence of slow particle transport in the many-body localized phase, Phys. Rev. B 102, 100202 (2020).
  * Šuntajs  _et al._ [2020] J. Šuntajs, J. Bonča, T. Prosen, and L. Vidmar, Ergodicity breaking transition in finite disordered spin chains, Phys. Rev. B 102, 064207 (2020).
  * Abanin  _et al._ [2021] 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, Annals of Physics 427, 168415 (2021).
  * Sels and Polkovnikov [2021] D. Sels and A. Polkovnikov, Dynamical obstruction to localization in a disordered spin chain, Phys. Rev. E 104, 054105 (2021).
  * Sels [2022] D. Sels, Bath-induced delocalization in interacting disordered spin chains, Phys. Rev. B 106, L020202 (2022).
  * Sierant and Zakrzewski [2022] P. Sierant and J. Zakrzewski, Challenges to observation of many-body localization, Phys. Rev. B 105, 224203 (2022).
  * Morningstar  _et al._ [2022] 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).
  * Crowley and Chandran [2022] 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).
  * Long  _et al._ [2023] D. M. Long, P. J. D. Crowley, V. Khemani, and A. Chandran, Phenomenology of the prethermal many-body localized regime, Phys. Rev. Lett. 131, 106301 (2023).
  * Evers  _et al._ [2023] F. Evers, I. Modak, and S. Bera, Internal clock of many-body delocalization, Phys. Rev. B 108, 134204 (2023).
  * Chávez  _et al._ [2023] D. A. Chávez, C. Artiaco, T. K. Kvorning, L. Herviou, and J. H. Bardarson, Ultraslow growth of number entropy in an l-bit model of many-body localization (2023), arXiv:2312.13420 [cond-mat.dis-nn] .
  * Zhang  _et al._ [2016] L. Zhang, V. Khemani, and D. A. Huse, A floquet model for the many-body localization transition, Phys. Rev. B 94, 224202 (2016).
  * Lezama  _et al._ [2019] T. L. M. Lezama, S. Bera, and J. H. Bardarson, Apparent slow dynamics in the ergodic phase of a driven many-body localized system without extensive conserved quantities, Phys. Rev. B 99, 161106 (2019).
  * Sierant  _et al._ [2023] P. Sierant, M. Lewenstein, A. Scardicchio, and J. Zakrzewski, Stability of many-body localization in floquet systems, Phys. Rev. B 107, 115132 (2023).
  * Kim and Huse [2013] H. Kim and D. A. Huse, Ballistic spreading of entanglement in a diffusive nonintegrable system, Phys. Rev. Lett. 111, 127205 (2013).
  * Burke  _et al._ [2023] P. C. Burke, G. Nakerst, and M. Haque, Assigning temperatures to eigenstates, Phys. Rev. E 107, 024102 (2023).
  * Lazarides  _et al._ [2014] A. Lazarides, A. Das, and R. Moessner, Periodic thermodynamics of isolated quantum systems, Phys. Rev. Lett. 112, 150401 (2014).
  * D’Alessio and Rigol [2014] L. D’Alessio and M. Rigol, Long-time behavior of isolated periodically driven interacting lattice systems, Phys. Rev. X 4, 041048 (2014).
  * Abanin  _et al._ [2015] D. A. Abanin, W. De Roeck, and F. m. c. Huveneers, Exponentially slow heating in periodically driven many-body systems, Phys. Rev. Lett. 115, 256803 (2015).
  * Mori  _et al._ [2016] T. Mori, T. Kuwahara, and K. Saito, Rigorous bound on energy absorption and generic relaxation in periodically driven quantum systems, Phys. Rev. Lett. 116, 120401 (2016).
  * Abanin  _et al._ [2017] D. A. Abanin, W. De Roeck, W. W. Ho, and F. m. c. Huveneers, Effective hamiltonians, prethermalization, and slow energy absorption in periodically driven many-body systems, Phys. Rev. B 95, 014112 (2017).
  * Mallayya and Rigol [2019] K. Mallayya and M. Rigol, Heating rates in periodically driven strongly interacting quantum many-body systems, Phys. Rev. Lett. 123, 240603 (2019).
  * Ponte  _et al._ [2015] P. Ponte, Z. Papić, F. m. c. Huveneers, and D. A. Abanin, Many-body localization in periodically driven systems, Phys. Rev. Lett. 114, 140401 (2015).
  * Lazarides  _et al._ [2015] A. Lazarides, A. Das, and R. Moessner, Fate of many-body localization under periodic driving, Phys. Rev. Lett. 115, 030402 (2015).
  * Rehn  _et al._ [2016] J. Rehn, A. Lazarides, F. Pollmann, and R. Moessner, How periodic driving heats a disordered quantum spin chain, Phys. Rev. B 94, 020201 (2016).
  * Bordia  _et al._ [2017] P. Bordia, H. Lüschen, U. Schneider, M. Knap, and I. Bloch, Periodically driving a many-body localized quantum system, Nature Physics 13, 460 (2017).
  * Haldar [2023] A. Haldar, Slow dynamics and kohlrausch relaxation in isolated disordered many-body systems (2023), arXiv:2302.12275 [cond-mat.stat-mech] .
  * Ikeda and Polkovnikov [2021] T. N. Ikeda and A. Polkovnikov, Fermi’s golden rule for heating in strongly driven floquet systems, Phys. Rev. B 104, 134308 (2021).
  * Rakcheev and Läuchli [2022] A. Rakcheev and A. M. Läuchli, Estimating heating times in periodically driven quantum many-body systems via avoided crossing spectroscopy, Phys. Rev. Res. 4, 043174 (2022).
  * Page [1993] D. N. Page, Average entropy of a subsystem, Phys. Rev. Lett. 71, 1291 (1993).
  * Lakshminarayan [2001] A. Lakshminarayan, Entangling power of quantized chaotic systems, Phys. Rev. E 64, 036207 (2001).
  * De Tomasi and Khaymovich [2020] G. De Tomasi and I. M. Khaymovich, Multifractality meets entanglement: Relation for nonergodic extended states, Phys. Rev. Lett. 124, 200602 (2020).
  * De Roeck and Huveneers [2017] W. De Roeck and F. m. c. Huveneers, Stability and instability towards delocalization in many-body localization systems, Phys. Rev. B 95, 155129 (2017).
  * Thiery  _et al._ [2018] T. Thiery, F. Huveneers, M. Müller, and W. De Roeck, Many-body delocalization as a quantum avalanche, Phys. Rev. Lett. 121, 140601 (2018).
  * Gopalakrishnan  _et al._ [2015] S. Gopalakrishnan, M. Müller, V. Khemani, M. Knap, E. Demler, and D. A. Huse, Low-frequency conductivity in many-body localized systems, Phys. Rev. B 92, 104202 (2015).
  * Garratt  _et al._ [2021] S. J. Garratt, S. Roy, and J. T. Chalker, Local resonances and parametric level dynamics in the many-body localized phase, Phys. Rev. B 104, 184203 (2021).
  * Garratt and Roy [2022] S. J. Garratt and S. Roy, Resonant energy scales and local observables in the many-body localized phase, Phys. Rev. B 106, 054309 (2022).
  * Léonard  _et al._ [2023] J. Léonard, S. Kim, M. Rispoli, A. Lukin, R. Schittko, J. Kwan, E. Demler, D. Sels, and M. Greiner, Probing the onset of quantum avalanches in a many-body localized system, Nature Physics 19, 481 (2023).
  * Bera  _et al._ [2017] 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).
  * Chanda  _et al._ [2020] 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).
  * Sierant  _et al._ [2020] P. Sierant, M. Lewenstein, and J. Zakrzewski, Polynomially filtered exact diagonalization approach to many-body localization, Phys. Rev. Lett. 125, 156601 (2020).
  * Doggen  _et al._ [2021] E. V. Doggen, I. V. Gornyi, A. D. Mirlin, and D. G. Polyakov, Many-body localization in large systems: Matrix-product-state approach, Ann. Phys. 435, 168437 (2021).

##  I Appendix

###  A Finite $L$ dependence

Figure 6 upper panel shows the system size dependence of the decay of the
disorder averaged inverse local temperature $\overline{\beta}(n)$. With
decreasing $\Gamma$, i.e., increasing the strength of the disorder, the
$L$-dependence is more severe for finite time simulations. The system
thermalizes quickly for larger $\Gamma=0.7$, and the data shows less finite-
size corrections.

###  B Estimation of decay time $\tau_{j}$

Figure 6 lower panel shows the decay of $\beta_{j}(n)$ for two typical
disorder configurations for different values of disorder strength, and $L=20$.
The dashed line represents the estimated curve using the decay time $\tau_{j}$
defined in Eq. (3). The estimate of $\tau_{j}$ reasonably reproduces the decay
for pure exponential traces. When $\beta_{j}(n)$ shows non-exponential decay
for a stronger disorder, the $\tau_{j}$ gives only a rough estimate as visible
at smaller $\Gamma\gtrsim 0.5$.

Figure 6: Upper panel: $L$-dependence of the decay of $\beta(n)$ with $n$ for
different values of $\Gamma=0.7,0.6,0.5$. Lower panel: Decay of $\beta_{j}(n)$
for two typical disorder configurations. The dashed lines show the curve
corresponding to the estimated decay times defined in Eq. (3).

###  C Proof of ${{\mathscr{S}}}_{\text{P}}=2{{\mathscr{S}}}_{\text{d}}$

Here, we establish the relation between the Fock space (FS) participation
entropy $\mathscr{S}_{\text{P}}$ and the diagonal entropy in the FS basis
defined in Eq. (4).

$\displaystyle\mathscr{S}_{\mathrm{P}}=-\sum_{i=1}^{2^{L}}|c_{i}(t)|^{2}\log(|c_{i}(t)|^{2})$

from the definition of the participation entropy. Now, the probabilities,
$|c_{i}(t)|^{2}$ can be written in terms of the diagonal elements of the
density matrix of the system, $\rho$ in the following way. We note,
$|c_{i}(t)|^{2}=c_{i}(t)c_{i}^{*}(t)=\langle
i|\psi\rangle\langle\psi|i\rangle=\langle i|\rho|i\rangle$, where
$\rho=|\psi\rangle\langle\psi|$ is the system’s density matrix.

$\mathscr{S}_{\mathrm{P}}=-\sum_{i=1}^{2^{L}}\rho_{ii}\log\rho_{ii}$

Here, we can decompose the i-th diagonal element of the full density matrix,
i.e., $\rho_{ii}$, as a product of diagonal elements of the sub-system density
matrices.

$\displaystyle\rho_{ii}$ $\displaystyle=\langle
i|\rho_{A}\otimes\rho_{B}|i\rangle$ $\displaystyle=\langle j|\langle
k|\rho_{A}\otimes\rho_{B}|j\rangle|k\rangle$
$\displaystyle=\rho_{jj}^{A}\rho_{kk}^{B}$ (5)

here, $|i\rangle=|j\rangle|k\rangle$. Therefore,

$\displaystyle\mathscr{S}_{\mathrm{P}}$
$\displaystyle=-\sum_{i=1}^{2^{L}}\rho_{ii}\log\rho_{ii}$
$\displaystyle=\sum_{j=1}^{2^{L/2}}\sum_{k=1}^{2^{L/2}}\rho_{jj}^{A}\rho_{kk}^{B}\log(\rho_{jj}^{A}\rho_{kk}^{B})\hskip
51.21495pt[\text{using Eq. \ref{decomposition of rho into rhoA and rhoB}}]$
$\displaystyle=\sum_{j=1}^{2^{L/2}}\sum_{k=1}^{2^{L/2}}\rho_{jj}^{A}\rho_{kk}^{B}\left(\log(\rho_{jj}^{A})+\log(\rho_{kk}^{B})\right)$
$\displaystyle=\sum_{k=1}^{2^{L/2}}\rho_{kk}^{B}\sum_{j=1}^{2^{L/2}}\rho_{jj}^{A}\log(\rho_{jj}^{A})+\sum_{j=1}^{2^{L/2}}\rho_{jj}^{A}\sum_{k=1}^{2^{L/2}}\rho_{kk}^{B}\log(\rho_{kk}^{B})$
$\displaystyle=2\sum_{j=1}^{2^{L/2}}\rho_{jj}^{A}\log(\rho_{jj}^{A})$
$\displaystyle=2\mathscr{S}_{\mathrm{d}}$ (6)