# Resource Allocation and Resolution Control in the Metaverse with Mobile
Augmented Reality

Peiyuan Si1, Jun Zhao1, Huimei Han1,2, Kwok-Yan Lam1, Yang Liu1,3
1School of Computer Science & Engineering, Nanyang Technological University,
Singapore
2College of Information Engineering, Zhejiang University of Technology,
Hangzhou, Zhejiang, P.R. China
3School of Information and Communication Engineering, Dalian University of
Technology, Dalian, China
<EMAIL_ADDRESS>{junzhao, huimei.han<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

With the development of blockchain and communication techniques, the Metaverse
is considered as a promising next-generation Internet paradigm, which enables
the connection between reality and the virtual world. The key to rendering a
virtual world is to provide users with immersive experiences and virtual
avatars, which is based on virtual reality (VR) technology and high data
transmission rate. However, current VR devices require intensive computation
and communication, and users suffer from high delay while using wireless VR
devices. To build the connection between reality and the virtual world with
current technologies, mobile augmented reality (MAR) is a feasible alternative
solution due to its cheaper communication and computation cost. This paper
proposes an MAR-based connection model for the Metaverse, and proposes a
communication resources allocation algorithm based on outer approximation (OA)
to achieve the best utility. Simulation results show that our proposed
algorithm is able to provide users with basic MAR services for the Metaverse,
and outperforms the benchmark greedy algorithm.

###### Index Terms:

Metaverse, resource allocation, mobile augmented reality, outer approximation

††publicationid: pubid:  978-1-6654-3540-6/22 © 2022 IEEE

##  I Introduction

The development of blockchain and communication techniques motivated intensive
interest in the Metaverse, which is considered a next-generation Internet
paradigm [1, 2, 3]. Attracted by the potential of the Metaverse, many
governments and companies around the world are planning and preparing for the
upcoming Metaverse era; e.g., South Korea is promoting lessons based on the
Metaverse, Facebook announced that it would become a Metaverse company and
renamed itself as Meta, and Tencent has invested in an AR platform called
“Avakin life” [4, 5, 6].

Current challenges. One of the key issues in the Metaverse is to connect the
virtual world and the real world with the support of extend reality (XR),
including virtual reality (VR) and augmented reality (AR) [7, 8, 9, 10].
Current research on the Metaverse mainly focuses on solving the problem of
communication and computation for VR to provide users with immersive
experience [11, 12]. However, mobile VR services require extremely high data
rate, which is not easy to achieve even under the context of 5G. Besides, VR
users suffer and have to bear with the high weight of current VR devices,
which is another problem to be solved.

Mobile augmented reality (MAR) is a possible alternative option for VR, and
also an important component of the Metaverse [13, 14, 15, 16]. Compared to VR,
MAR combines reality and the virtual world rather than creating a fully
virtual world, which saves communication and computational cost [17, 18, 19].
Besides, current AR devices have a significant advantage over VR in weight,
making them more comfortable and safer as wearable devices. AR also has its
unique advantage in some applications, e.g., navigation, health care, tourism,
shopping and education, where the interaction with reality is required [20,
21, 22].

Related work. Although AR requires lower data rate than VR, efficient
allocation of communication resources is still necessary due to the massive
number of users and devices connected to the Metaverse server. To improve the
communication resource efficiency and quality of service (QoS), MEC and
reinforcement learning (RL) for VR/AR service have attracted much attention
[23, 24, 25, 26, 27, 28, 29]. Feng _et al._ [28] proposed a smart VR
transmission mode scheme based on RL to optimize the D2D system throughput and
achieve a balance between performance and resource efficiency. Chen _et al._
[23] introduced an RL-based energy-efficient MEC framework for AR services
with task offloading and resource allocation to release the burden at the
terminal. A recent work [30] applies deep RL to MAR services of the metaverse
over 6G wireless networks. Resolution control is also one of the solutions to
improve resource efficiency [31]. Higher resolution brings better QoS at the
cost of occupying more communication resources while lower resolution improves
resource efficiency. Thus, finding the balance between QoS and communication
resources such as power and bandwidth is of vital importance to Metaverse MAR
service.

In this paper, we propose an MAR-based Metaverse model with resolution control
and resource allocation algorithm. In our proposed resource allocation
optimization problem, both QoS and power consumption are included in the
utility function for the balance between user experience and the energy
efficiency of the service provider. To solve the mixed-integer nonlinear
programming (MINLP) optimization problem, we propose a resource allocation
algorithm based on outer approximation (OA), which guarantees global optimum.
Simulation results show that our algorithm is able to maximize the utility
with given communication resources, and outperforms the benchmark greedy
algorithm.

Contributions. The contributions of this paper are as follows:

  * •

A MAR-based Metaverse model is proposed for applications where interaction
with reality is required.

  * •

A resolution control and resource allocation optimization problem for
Metaverse MAR services is proposed and solved by an OA algorithm.

  * •

Simulation of the proposed algorithm is implemented with the comparison to a
benchmark greedy algorithm under various parameter settings. The results show
the advantage of our proposed algorithm and the feasibility of the proposed
Metaverse MAR service model.

The rest of this paper is organized as follows. Section II introduces the
proposed system model. The problem formulation and solution are presented in
Section III and Section IV respectively. Section V shows the simulation
results and the corresponding explanation. The conclusion of this paper is
discussed in Section VI.

##  II System Model and Problem Formulation

Figure 1: System model of the Metaverse with mobile augmented reality.

In our proposed Metaverse MAR service model, we consider MAR service for
multiple users supported by a base station in a particular area, which is
connected to Metaverse servers. As shown in Fig. 1, MAR users can get
Metaverse services with data synchronized to their avatars in the Metaverse
while interacting with the real world with the support of the base station,
which controls the MAR service quality by switching the resolution, e.g.,
480P, 720P, 1080P. Due to the limit of total transmission power, the base
station needs to dynamically allocate power for users and control the
resolution of MAR service according to the channel condition to achieve the
best utility.

We assmue that the base station provides $N$ users with MAR service through
$N$ channels. With fixed bandwidth and limited transmit power, the base
station allocates transmit power to $N$ users, and the transmission rate for
user $n$ is given by

$\displaystyle{{r}_{n}}=B{{\log}_{2}}\left(1+\frac{{{p}_{n}}{{g}_{n}}}{{{\sigma}^{2}}}\right),$
(1)

where $B$ and $p_{n}$ denotes the bandwidth and transmit power for user $n$
respectively, and $\sigma^{2}$ denotes the power of Gaussian noise. $g_{n}$
denotes the channel gain between base station and user $n$. Let
$\boldsymbol{p}:=[p_{1},p_{2},\ldots,p_{N}]$. The base station is able to
provide users with different resolutions of MAR service according to the
channel condition. The required transmission rate for real time MAR service is
denoted by

$\displaystyle{{C}_{n}}\in\\{{{C}_{1}},{{C}_{2}},{{C}_{3}}\\},$ (2)

where ${{C}_{1}},{{C}_{2}},{{C}_{3}}$ denote the required transmission rate of
360P, 720P and 1080P MAR service respectively. The base station can provide
each user with only one type of service, and the selection indicator is given
by

$\displaystyle\boldsymbol{I}=\\{{{I}_{n,1}},{{I}_{n,2}},{{I}_{n,3}}|n\in[N]\\},$
(3)

where ${{I}_{n,1}},{{I}_{n,2}},{{I}_{n,3}}\in\\{0,1\\}$,
$\sum_{i=1}^{3}{{I}_{n,i}}=1$ and $[N]:=\\{1,2,..,N\\}$. To evaluate the
quality of service (QoS) for users with given service selection, we introduce
a QoS factor $Q_{n}$, which is given by

$\displaystyle{{Q}_{n}}={{\left(\frac{{{C}_{n}}}{{{R}_{th}}}\right)}^{\gamma}}=\sum\limits_{i=1}^{3}{{{I}_{n,i}}{{\left(\frac{{{C}_{i}}}{{{R}_{th}}}\right)}^{\gamma}}},$
(4)

where $R_{th}$ denotes the constant reference service quality, and $\gamma$
denotes the exponential factor of QoS. The best QoS is given by

$\displaystyle{{Q}_{0}}={{\left(\frac{{{C}_{3}}}{{{R}_{th}}}\right)}^{\gamma}},$
(5)

QoS only evaluates the experience of individual user, which does not fully
represent the utility of the resource allocation scheme. If the base station
provides more users with higher QoS under given bandwidth resource, more power
will be occupied to reach the required transmission rate. Thus, we take the
total power consumption into consideration for the utility function, which is
given by

$\displaystyle
U=(1-\lambda-\mu){{\eta}_{q}}\sum\limits_{n=1}^{N}{{{Q}_{n}}}-\lambda{{\eta}_{p}}\sum\limits_{n=1}^{N}{{{p}_{n}}}+\mu{{\eta}_{r}}\sum\limits_{n=1}^{N}{\Delta{{r}_{n}}},$
(6)
$\displaystyle\Delta{{r}_{n}}=\sum\limits_{i=1}^{3}{{{C}_{i}}I_{n,i}^{(0)}}-B{{\log}_{2}}\left(1+\frac{{{p}_{n}}{{g}_{n}}}{{{\sigma}^{2}}}\right),$
(7)

where ${{\eta}_{q}}={{\left(N{{Q}_{0}}\right)}^{-1}}$,
${{\eta}_{r}}={{{{C}_{2}}}^{-1}}$ and ${{\eta}_{p}}={{P}^{-1}}$ are
normalization factors, $\lambda$, $\mu$ and $1-\lambda-\mu$ denote the levels
of concern for energy consumption, redundancy of transmission rate and QoS
respectively. We consider the redundancy of transmission rate as a positive
factor because it provides robustness against the turbulence of channel. To
maximize the utility of the resource allocation scheme under given total
transmission power limit, the optimization problem is given by

$\displaystyle
P1:\underset{\boldsymbol{I},\boldsymbol{p}}{\mathop{\min}}\,\text{
}\lambda{{\eta}_{p}}\sum\limits_{n=1}^{N}{{{p}_{n}}-\mu{{\eta}_{r}}\sum\limits_{n=1}^{N}{\Delta{{r}_{n}}}-}(1-\lambda-\mu){{\eta}_{q}}\sum\limits_{n=1}^{N}{{{Q}_{n}}}$
$\displaystyle subject\text{ }to$ $\displaystyle
Con1:\sum\limits_{i=1}^{3}{{{C}_{i}}I_{n,i}^{(0)}}-B{{\log}_{2}}\left(1+\frac{{{p}_{n}}{{g}_{n}}}{{{\sigma}^{2}}}\right)\leq
0,\forall n\in[N],$ $\displaystyle Con2:\sum\limits_{n=1}^{N}{{{p}_{n}}}\leq
P,$ $\displaystyle Con3:\sum\limits_{i=1}^{3}{I_{n,i}}=1,\forall
n\in[N],I_{n,i}\in\\{0,1\\},$ (8)

where constraint $Con1$ indicates that data transmission rate should meet the
requirement of the selected service quality. Constraint $Con2$ denotes the
total power constraint, and constraint $Con3$ indicates that the base station
can provide each user with only one type of service. Problem $P1$ is a mixed-
integer nonlinear programming (MINLP) problem, which can be solved by outer
approximating (OA) method. The solution of $P1$ is presented in Section IV.

##  III Problem Solution

The MINLP problem $P1$ can be solved by OA method in different steps, which
are explained in the subsections below.

###  III-A Solve the Nonlinear Programming Problem with Given Integer
Variables

At the beginning of the OA algorithm, i.e., the first iteration, we give
initial values to the integer variables matrix $\boldsymbol{I}$ in the
feasible region, which are denoted as $\boldsymbol{I}^{(0)}$. Then substitute
$\boldsymbol{I}^{(0)}$ into $P1$ to formulate the nonlinear programming
problem $P1.1$, which is given by

$\displaystyle P1.1:\underset{\boldsymbol{p}}{\mathop{\min}}\,\text{
}\lambda{{\eta}_{p}}\sum\limits_{n=1}^{N}{{{p}_{n}}-\mu{{\eta}_{r}}\sum\limits_{n=1}^{N}{\Delta{{r}_{n}}}}-{{A}_{1}}$
$\displaystyle subject\text{ }to$ $\displaystyle
Con4:B{{\log}_{2}}\left(1+\frac{{{p}_{n}}{{g}_{n}}}{{{\sigma}^{2}}}\right)-\sum\limits_{i=1}^{3}{{{C}_{i}}I_{n,i}^{(0)}}\geq
0,\forall n\in[N],$ $\displaystyle Con5:P-\sum\limits_{n=1}^{N}{{{p}_{n}}}\geq
0,$ (9)

where the constant $A_{1}$ is given by

$\displaystyle{{A}_{1}}=\text{
}(1-\lambda-\mu){{\eta}_{q}}\sum\limits_{n=1}^{N}{\sum\limits_{i=1}^{3}{I_{n,i}^{(0)}{{\left(\frac{{{C}_{i}}}{{{R}_{th}}}\right)}^{2}}}}.$
(10)

With given integer variables, problem $P1$ is transformed into nonlinear
programming (NLP) $P1.1$, which can be solved by interior point method, whose
worst case iteration complexity in this paper is
${\mathrm{O}}\left(\sqrt{n}\log n\log\frac{n}{\epsilon}\right)$ [32].

We introduce an auxiliary variable $s_{n}$ to reformulate the problem as

$\displaystyle P1.2:\underset{\boldsymbol{p}}{\mathop{\min}}\,\text{
}\lambda{{\eta}_{p}}\sum\limits_{n=1}^{N}{{{p}_{n}}-\mu{{\eta}_{r}}\sum\limits_{n=1}^{N}{\Delta{{r}_{n}}}}-{{A}_{1}}$
$\displaystyle subject\text{ }to$ $\displaystyle
Con4:B{{\log}_{2}}\left(1+\frac{{{p}_{n}}{{g}_{n}}}{{{\sigma}^{2}}}\right)-\sum\limits_{i=1}^{3}{{{C}_{i}}I_{n,i}^{(0)}}={{s}_{n}},\forall
n\in[N],$ $\displaystyle Con5:P-\sum\limits_{n=1}^{N}{{{p}_{n}}}\geq 0,$
$\displaystyle Con6:{{s}_{n}}\geq 0,\forall n\in[N].$ (11)

The barrier function is given by

$\displaystyle J(p,\tau)$
$\displaystyle=\lambda{{\eta}_{p}}\sum\limits_{n=1}^{N}{{{p}_{n}}-\mu{{\eta}_{r}}\sum\limits_{n=1}^{N}{\Delta{{r}_{n}}}}-{{A}_{1}}$
$\displaystyle\quad-\tau\left(\log\left({{s}_{n}}\right)+\log\left(P-\sum\limits_{n=1}^{N}{{{p}_{n}}}\right)\right),$
(12)

where $\tau$ is a positive barrier parameter. The perturbed KKT conditions for
$n\in[N]$ are given by [33]

$\displaystyle\nabla J(p,\tau)-\Delta{{r}_{n}}{{z}_{n}}=0,$ (13)
$\displaystyle\Delta{{r}_{n}}-{{s}_{n}}=0,$ (14) $\displaystyle
s_{n}z_{n}=\tau,$ (15) $\displaystyle s_{n}\geq 0,z_{n}\geq 0,$ (16)

where $z_{n},n\in[N]$ are Lagrange multiplier-inspired dual variables. The
interior point method starts with a feasible point which satisfies the
perturbed KKT conditions with small $\tau$, and continue to find smaller
$\tau$. As the value of $\tau$ approaches zero, the solution is expected to
converge to a point which satisfies the KKT point of problem $P1.1$.

With the obtained optimal power allocation, the solution of $P1.1$ at $t$-th
iteration is denoted as $z_{ub}^{(t)}\left(\boldsymbol{p}_{ub}^{(t)}\right)$,
which is a lower bound of the global optimum, where
$\boldsymbol{p}_{ub}^{(t)}$ denotes the optimal solution at the $t$-th
iteration.

###  III-B Solve the Mixed-integer Linear Programming Problem

After obtaining the optimal solution for $P1.2$, we formulate the first-order
Taylor expansion of the original problem $P1$ at point
$\boldsymbol{p}_{ub}^{(t)}$ as

$\displaystyle P1.3:\underset{\boldsymbol{I,p}}{\mathop{\min}}\,\text{
}\lambda{{\eta}_{p}}\sum\limits_{n=1}^{N}{{{p}_{n}}}-\mu{{\eta}_{r}}\sum\limits_{n=1}^{N}{\left({{A}_{2}}-\sum\limits_{i=1}^{3}{{{C}_{i}}{{I}_{n,i}}}\right)}$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-(1-\lambda-\mu){{\eta}_{q}}\sum\limits_{n=1}^{N}{{{Q}_{n}}}$
$\displaystyle subject\text{ }to$ $\displaystyle
Con7:{{A}_{2}}-\sum\limits_{i=1}^{3}{{{C}_{i}}{{I}_{n,i}}}\leq 0,\forall
n\in[N],$ $\displaystyle Con2,Con3,$ (17)

where

$\displaystyle{{A}_{2}}=B{{\log}_{2}}\left(1+\frac{p_{n}^{(t)}{{g}_{n}}}{{{\sigma}^{2}}}\right)+\frac{B{{g}_{n}}}{\ln
2\left({{\sigma}^{2}}+p_{n}^{(t)}{{g}_{n}}\right)}\left({{p}_{n}}-p_{n}^{(t)}\right).$
(18)

Problem $P1.3$ is a mixed-integer linear programming program (MILP), which can
be solved by existing MILP solvers with branch-and-bound algorithm. As MILP
problems are NP-hard, there are no polynomial complexity algorithms, and thus
the complexity in the worst case is still exponential. However, the branch-
and-bound algorithm reduces the average complexity compared to brute force
search [34]. The precise estimation of branch-and-bound algorithm complexity
requires the probability that a node in the branch-and-bound tree is optimal,
which is hard to obtain in the optimization problem of this paper. To analyze
the complexity of this algorithm, we implemented running time experiments
which show acceptable results that the algorithm is able to converge within 30
seconds with CPU i7-9750H and 16GB memory.

The optimal solution of problem $P1.3$ at the $t$-th iteration is denoted as
$z_{lb}^{(t)}\left(\boldsymbol{I}_{lb}^{(t)},\boldsymbol{p}_{lb}^{(t)}\right)$,
which is a lower bound of the global optimum.

###  III-C Compare the Gap Between $z_{lb}^{(t)}$ and $z_{ub}^{(t)}$

As the lower bound and upper bound of global optimum are obtained, the gap
between them is given by

$\displaystyle G=\left|z_{ub}^{(t)}-z_{lb}^{(t)}\right|.$ (19)

The algorithm is considered to be converged when $G\leq\epsilon$, where
$\epsilon$ is a given precision factor. If $G\geq\epsilon$, substitute the
integer variables $\boldsymbol{I}_{lb}^{(t)}$ into the original problem $P1$
for the next iteration. The overall algorithm is shown in Algorithm 1.

Algorithm 1 Outer Approximating Algorithm for Communication Resource
Allocation in the Metaverse

Initialize: $t=1$
Input: ${\textbf{I}^{(0)}}$
Step 1.1: Substitute ${\boldsymbol{I}^{(t-1)}}$ into $P1$ to formulate
nonlinear programming (NLP) problem $P1.1$
Step 1.2: Solve $P1.1$ by interior method to get lower bound
$z_{ub}^{(t)}\left(\boldsymbol{p}_{ub}^{(t)}\right)$ of the global optimum
Step 2.1: Formulate the MILP problem $P1.3$ with the first-order Taylor
expansion at point $\boldsymbol{p}_{ub}^{(t)}$
Step 2.2: Solve $P1.3$ by MILP solvers to get
$z_{lb}^{(t)}\left(I_{lb}^{(t)},p_{lb}^{(t)}\right)$
Step 3: If $G=\left|z_{ub}^{(t)}-z_{lb}^{(t)}\right|\leq\epsilon$, Output:
$\\{{\boldsymbol{I}^{(t)}},{{\boldsymbol{p}}^{(t)}}\\}$, Else Go to Step 1.1

##  IV Simulation Results

In this section, we present the simulation results of the proposed OA-based
resource allocation algorithm and the comparison with a benchmark greedy
algorithm. In the simulation, the setting of constant parameters is given in
Table I, and the benchmark greedy algorithm is given in Algorithm 2. The
channel gain in simulation is set as

$\displaystyle{{g}_{n}}$
$\displaystyle=20{{\log}_{10}}({{d}_{n}})+20{{\log}_{10}}(f)-147.55,$ (20)
$\displaystyle{{d}_{n}}$ $\displaystyle={{D}_{f}}{{d}_{0}}(n),$ (21)

where ${{d}_{0}}(n)$ denotes the reference distance between base station and
user $n$, and ${{D}_{f}}$ denotes the distance multiply factor. The default
value of $\mu$ is set to 0.1.

TABLE I: Constant Parameter Setting Parameter and Physical Meaning | Value
---|---
Exponential factor of QoS ($\gamma$) | $2$
Number of users ($N$) | $5$
Precision factor ($\epsilon$) | $10^{-3}$
Bandwidth ($B$) |  $5$MHz
Frequency ($f$) |  $28$GHz (i.e., 5G spectrum)
Reference service quality ($R_{th}$) |  $0.77$Mbps
Power of Gaussian noise ($\sigma^{2}$) |  $5\times{{10}^{-8}}$W
Reference distance ($d_{0}$) |  $\left(5,4,3,2,1\right)$(m)
Required transmission rates $\left(C_{1},C_{2},C_{3}\right)$ |  $\left(0.77,1.92,3.84\right)$(Mbps)
Algorithm 2 Benchmark Greedy Algorithm

Initialize: $P_{rest}=P$, $U_{max}=-\infty$
For $n\leq N$
For $i\leq 3$
$I_{n,i}^{(t)}=1$, $p_{n}^{(t)}$ = minimum required power
Calculate the utility $U$ and remaining power $P_{rest}$ under current
selection
If $U>U_{max}$ and $P_{rest}\geq 0$
Update $U_{max}=U$, $I_{n,i}=I_{n,i}^{(t)}$, $p_{n}=p_{n}^{(t)}$
Output: $U_{max}$, $I_{n,i}$, $p_{n}$

Figure 2: Utility comparison between two algorithms with different $D_{f}$ and
$\lambda$.

Fig. 2 shows the utility comparison between two algorithms with different
$D_{f}$ and $\lambda$. We can find that the utility of both algorithms
decrease with increasing user-to-BS distance because the base station needs
more power to provide the same service under smaller channel power gain. The
OA algorithm outperforms the greedy algorithm in two different cases, i.e.,
$\lambda=0.1$ and $\lambda=0.5$, because it ensures global optimum while the
greedy algorithm only obtains a local optimum.

Figure 3: Utility comparison between two algorithms with different $P$ and
$\lambda$.

Fig. 3 shows the utility comparison between two algorithms with different $P$
and $\lambda$. The utilities of both algorithms increase as the total
available transmission power $P$ increases because less power is required for
the same service, and the base station is able to provide users with better
services and higher QoS. Fig. 3 also shows the greedy algorithm’s ability to
find the local optimum and the advantage of OA algorithm in more general cases
due to the guarantee of the global optimum.

Figure 4: Power allocation with $\lambda=\mu=0.1$, $D_{f}=5$, $P=50$W. Figure
5: Service selection with $\lambda=\mu=0.1$, $D_{f}=5$, $P=50$W.

Fig. 4 and Fig. 5 show the power allocation and service selection under given
$D_{f}=5$, $\lambda=0.1$ and $P=50$W. The powers allocated by the greedy
algorithm to users 1, 2 and 3 decrease as the index increase because the
distances between the base station and users are $d_{3}<d_{2}<d_{1}$. As the
distance decreases, the base station needs less power for basic 360P service.
Due to the relatively poor channel condition of users 1, 2 and 3, it is
expensive to improve the QoS of these three users. The OA algorithm results in
a different policy for power allocation due to its consideration of robustness
against channel turbulence, which is also shown in 5. Although user 1 and user
2 are assigned with same service quality, user 1 requires more power due to
the longer distance to the base station, and it is the same case for user 3
and user 4. User 5 is assigned with large power because the robustness against
channel turbulence contributes to the utility function.

Figure 6: Convergence of OA algorithm.

Fig. 6 shows the convergence of the proposed OA algorithm for resource
allocation under two different sets of parameters. We can find that the upper
bound $z_{ub}^{(t)}$ and lower bound $z_{lb}^{(t)}$ converge within 5
iterations regardless of the parameter setting. In the simulation, we set the
convergence threshold $\epsilon=10^{-3}$ to limit the gap between
$z_{ub}^{(t)}$ and $z_{lb}^{(t)}$, but in most cases the gap becomes zero
after several iterations, which further guarantees the convergence.

Figure 7: Influence of $\gamma$ on service selection.

Fig. 7 shows the influence of $\gamma$ on service selection under $D_{f}=5$,
$P=50$W, $\lambda=0.1$, where the base station takes the power consumption as
the main concern, and the power resource is adequate to provide all users with
the best service. In this case, the exponential factor of QoS has a great
influence on the policy decisions of the base station. With larger $\gamma$,
the base station gets more benefit in utility through improving the QoS, and
becomes more willing to provide the users with better service even if it gets
more penalty from energy consumption. Thus, setting $\gamma$ properly is
critical to the balance between power consumption and QoS. As $\gamma$
increases, the base station tends to improve QoS of the users with larger
index first because they are closer to the base station than those with
smaller index, which requires less power to achieve the same QoS.

##  V Conclusion

In this paper, we proposed a MAR-based Metaverse model for applications that
require interaction with the real world. In order to maximize the utility of
the base station, we formulate a resolution control and resource allocation
optimization problem, and solve the MINLP problem with OA algorithm. The
simulation results indicate that our proposed resource allocation algorithm
outperforms the benchmark greedy algorithm due the guarantee of global
optimum. Our work also shows the feasibility of MAR to provide basic Metaverse
service and the advantage of our proposed algorithm over benchmark greedy
algorithm.

## Acknowledgement

This research/project is supported by the National Research Foundation,
Singapore under its Strategic Capability Research Centres Funding Initiative.
Any opinions, findings and conclusions or recommendations expressed in this
material are those of the author(s) and do not reflect the views of National
Research Foundation, Singapore. This research/project is supported by the
Ministry of Education, Singapore, under its Grant Tier 1 RG97/20, and Grant
Tier 1 RG24/20; in part by the NTU-Wallenberg AI, Autonomous Systems and
Software Program (WASP) Joint Project; and in part by the Singapore NRF
National Satellite of Excellence, Design Science and Technology for Secure
Critical Infrastructure under Grant NSoE DeST-SCI2019-0012. This work was also
supported in part by National Natural Science Foundation of China under Grants
62001419 and 62131016, in part by Zhejiang Provincial Natural Science
Foundation of China under Grant LQ21F010012. This work is also funded by
Dalian University of Technology, China under Grant No. DUT20RC(3)029.

## References

  * [1] M. Xu, D. Niyato, J. Kang, Z. Xiong, C. Miao, and D. I. Kim, “Wireless edge-empowered metaverse: A learning-based incentive mechanism for virtual reality,” _arXiv preprint arXiv:2111.03776_ , 2021.
  * [2] Y. Han, D. Niyato, C. Leung, C. Miao, and D. I. Kim, “A dynamic resource allocation framework for synchronizing metaverse with iot service and data,” _arXiv preprint arXiv:2111.00431_ , 2021.
  * [3] R. Cheng, N. Wu, S. Chen, and B. Han, “Will metaverse be nextg internet? vision, hype, and reality,” _arXiv preprint arXiv:2201.12894_ , 2022.
  * [4] H. Ning, H. Wang, Y. Lin, W. Wang, S. Dhelim, F. Farha, J. Ding, and M. Daneshmand, “A survey on metaverse: the state-of-the-art, technologies, applications, and challenges,” _arXiv preprint arXiv:2111.09673_ , 2021.
  * [5] S. M. Park and Y. G. Kim, “A metaverse: taxonomy, components, applications, and open challenges,” _IEEE Access_ , 2022.
  * [6] B. Egliston and M. Carter, “Critical questions for facebook’s virtual reality: Data, power and the metaverse,” _Internet Policy Review_ , vol. 10, no. 4, pp. 1–23, 2021.
  * [7] L. H. Lee, T. Braud, P. Zhou, L. Wang, D. Xu, Z. Lin, A. Kumar, C. Bermejo, and P. Hui, “All one needs to know about metaverse: A complete survey on technological singularity, virtual ecosystem, and research agenda,” _arXiv preprint arXiv:2110.05352_ , 2021.
  * [8] F. Hu, Y. Deng, W. Saad, M. Bennis, and A. H. Aghvami, “Cellular-connected wireless virtual reality: Requirements, challenges, and solutions,” _IEEE Communications Magazine_ , vol. 58, no. 5, pp. 105–111, 2020.
  * [9] I. Akyildiz and H. Guo, “Wireless extended reality (xr): Challenges and new research directions,” in _Journal Future and Evolving Technologies_. ITU, 2022.
  * [10] C. Chaccour, M. N. Soorki, W. Saad, M. Bennis, and P. Popovski, “Can terahertz provide high-rate reliable low latency communications for wireless vr?” _IEEE Internet of Things Journal_ , 2022.
  * [11] Y. Zhou, C. Pan, P. L. Yeoh, K. Wang, M. Elkashlan, B. Vucetic, and Y. Li, “Communication-and-computing latency minimization for uav-enabled virtual reality delivery systems,” _IEEE Transactions on Communications_ , vol. 69, no. 3, pp. 1723–1735, 2020.
  * [12] X. Liu, X. Li, and Y. Deng, “Learning-based prediction and proactive uplink retransmission for wireless virtual reality network,” _IEEE Transactions on Vehicular Technology_ , vol. 70, no. 10, pp. 10 723–10 734, 2021\.
  * [13] Z. Huang and V. Friderikos, “Proactive edge cloud optimization for mobile augmented reality applications,” in _2021 IEEE Wireless Communications and Networking Conference (WCNC)_. IEEE, 2021, pp. 1–6.
  * [14] C. Bermejo and P. Hui, “A survey on haptic technologies for mobile augmented reality,” _ACM Computing Surveys (CSUR)_ , vol. 54, no. 9, pp. 1–35, 2021\.
  * [15] Q. Liu, S. Huang, J. Opadere, and T. Han, “An edge network orchestrator for mobile augmented reality,” in _IEEE INFOCOM 2018-IEEE Conference on Computer Communications_. IEEE, 2018, pp. 756–764.
  * [16] J. Grubert, T. Langlotz, S. Zollmann, and H. Regenbrecht, “Towards pervasive augmented reality: Context-awareness in augmented reality,” _IEEE transactions on visualization and computer graphics_ , vol. 23, no. 6, pp. 1706–1724, 2017.
  * [17] Y. Siriwardhana, P. Porambage, M. Liyanage, and M. Ylianttila, “A survey on mobile augmented reality with 5g mobile edge computing: architectures, applications, and technical aspects,” _IEEE Communications Surveys & Tutorials_, vol. 23, no. 2, pp. 1160–1192, 2021.
  * [18] M. Erol-Kantarci and S. Sukhmani, “Caching and computing at the edge for mobile augmented reality and virtual reality (ar/vr) in 5g,” _Ad Hoc Networks_ , pp. 169–177, 2018.
  * [19] J. Ren, Y. He, G. Huang, G. Yu, Y. Cai, and Z. Zhang, “An edge-computing based architecture for mobile augmented reality,” _IEEE Network_ , vol. 33, no. 4, pp. 162–169, 2019.
  * [20] N. Xi, J. Chen, F. Gama, M. Riar, and J. Hamari, “The challenges of entering the metaverse: An experiment on the effect of extended reality on workload,” _Information Systems Frontiers_ , pp. 1–22, 2022.
  * [21] D. Buhalis and N. Karatay, “Mixed reality (mr) for generation z in cultural heritage tourism towards metaverse,” in _ENTER22 e-Tourism Conference_. Springer, 2022, pp. 16–27.
  * [22] C. Moro, C. Phelps, P. Redmond, and Z. Stromberga, “Hololens and mobile augmented reality in medical and health science education: A randomised controlled trial,” _British Journal of Educational Technology_ , vol. 52, no. 2, pp. 680–694, 2021.
  * [23] X. Chen and G. Liu, “Energy-efficient task offloading and resource allocation via deep reinforcement learning for augmented reality in mobile edge networks,” _IEEE Internet of Things Journal_ , vol. 8, no. 13, pp. 10 843–10 856, 2021.
  * [24] Y. Sun, Z. Chen, M. Tao, and H. Liu, “Communications, caching, and computing for mobile virtual reality: Modeling and tradeoff,” _IEEE Transactions on Communications_ , vol. 67, no. 11, pp. 7573–7586, 2019.
  * [25] Q. Liu and T. Han, “Dare: Dynamic adaptive mobile augmented reality with edge computing,” in _2018 IEEE 26th International Conference on Network Protocols (ICNP)_. IEEE, 2018, pp. 1–11.
  * [26] F. Guo, F. R. Yu, H. Zhang, H. Ji, V. C. Leung, and X. Li, “An adaptive wireless virtual reality framework in future wireless networks: A distributed learning approach,” _IEEE Transactions on Vehicular Technology_ , vol. 69, no. 8, pp. 8514–8528, 2020.
  * [27] C. Wang, S. Zhang, Y. Chen, Z. Qian, J. Wu, and M. Xiao, “Joint configuration adaptation and bandwidth allocation for edge-based real-time video analytics,” in _IEEE INFOCOM 2020-IEEE Conference on Computer Communications_. IEEE, 2020, pp. 257–266.
  * [28] L. Feng, Z. Yang, Y. Yang, X. Que, and K. Zhang, “Smart mode selection using online reinforcement learning for vr broadband broadcasting in d2d assisted 5g hetnets,” _IEEE Transactions on Broadcasting_ , vol. 66, no. 2, pp. 600–611, 2020.
  * [29] X. Liu and Y. Deng, “Learning-based prediction, rendering and association optimization for mec-enabled wireless virtual reality (vr) networks,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 10, pp. 6356–6370, 2021.
  * [30] T. J. Chua, W. Yu, and J. Zhao, “Resource allocation for mobile metaverse with the Internet of Vehicles over 6G wireless communications: A deep reinforcement learning approach,” in _8th IEEE World Forum on the Internet of Things (WFIoT)_ , 2022.
  * [31] J. Ahn, J. Lee, S. Yoon, and J. K. Choi, “A novel resolution and power control scheme for energy-efficient mobile augmented reality applications in mobile edge computing,” _IEEE Wireless Communications Letters_ , vol. 9, no. 6, pp. 750–754, 2019.
  * [32] M. R. Peyghami and S. F. Hafshejani, “Complexity analysis of an interior-point algorithm for linear optimization based on a new proximity function,” _Numerical Algorithms_ , vol. 67, no. 1, pp. 33–48, 2014.
  * [33] M. Mehrjoo, S. Moazeni, and X. Shen, “Resource allocation in ofdma networks based on interior point methods,” _Wireless Communications and Mobile Computing_ , vol. 10, no. 11, pp. 1493–1508, 2010.
  * [34] N. Thakoor, V. Devarajan, and J. Gao, “Computation complexity of branch-and-bound model selection,” in _2009 IEEE 12th International Conference on Computer Vision_. IEEE, 2009, pp. 1895–1900.


# Energy-recurrence Breakdown and Chaos in Disordered Fermi-Pasta-Ulam-Tsingou
Lattices

Zulkarnain  H. Susanto<EMAIL_ADDRESS>C. G. Antonopoulos  Department
of Mathematical Sciences, University of Essex, Wivenhoe Park, Colchester CO4
3SQ, United Kingdom  Department of Mathematics, Faculty of Mathematics and
Natural Sciences, Universitas Riau, Pekanbaru 28293, Indonesia  Department of
Mathematics, Khalifa University, PO Box 127788, Abu Dhabi, United Arab
Emirates

###### Abstract

In this paper, we consider the classic Fermi-Pasta-Ulam-Tsingou system as a
model of interacting particles connected by harmonic springs with a quadratic
nonlinear term (first system) and a set of second-order ordinary differential
equations with variability (second system) that resembles Hamilton’s equations
of motion of the Fermi-Pasta-Ulam-Tsingou system. In the absence of
variability, the second system becomes Hamilton’s equations of motion of the
Fermi-Pasta-Ulam-Tsingou system (first system). Variability is introduced to
Hamilton’s equations of motion of the Fermi-Pasta-Ulam-Tsingou system to take
into account inherent variations (for example, due to manufacturing
processes), giving rise to heterogeneity in its parameters. We demonstrate
that a percentage of variability smaller than a threshold can break the well-
known energy recurrence phenomenon and induce localization in the energy
normal-mode space. However, percentage of variability larger than the
threshold may make the trajectories of the second system blow up in finite
time. Using a multiple-scale expansion, we derive analytically a two normal-
mode approximation that explains the mechanism for energy localization and
blow up in the second system. We also investigate the chaotic behavior of the
two systems as the percentage of variability is increased, utilising the
maximum Lyapunov exponent and Smaller Alignment Index. Our analysis shows that
when there is almost energy localization in the second system, it is more
probable to observe chaos, as the number of particles increases.

###### keywords:

Fermi-Pasta-Ulam-Tsingou (FPUT) Hamiltonian , Chaos, Blow up , Maximum
Lyapunov exponent, Smaller Alignment Index (SALI), Multiple-scale expansion,
Two normal-mode approximation , Bifurcation analysis

###### MSC:

[2010] 00-01, 99-00

††journal:

##  1 Introduction

Debye suggested that thermal conductivity in a crystal is a consequence of
atom vibrations in the lattice [1, 2]. To model thermalization processes in
physical media, Fermi, Pasta, Ulam, with Tsingou’s help running the computer
simulations [3], considered a system of particles connected by harmonic
springs with a quadratic nonlinear term, i.e., the so-called FPUT lattice,
that is fixed at both ends.

A purely linear dynamics of the springs keeps energy, given to a single normal
mode, localized in that mode. However, introducing nonlinear interactions, one
would expect that energy introduced to one normal mode, would slowly spread to
other normal modes, until the system reaches a state of equipartition of
energy, i.e., the system relaxes to a thermal equilibrium. Contrary to this
expectation, Fermi, Pasta, Ulam, and Tsingou observed in their seminal paper
[3] in 1955 recurrences of the energy to its initial state, the so-called FPUT
recurrences, a phenomenon that led to numerous discoveries in mathematics and
physics [4, 5, 6, 7, 8] thereafter.

The recurrence phenomenon was explained by Zabusky and Kruskal [9] in real
space, who derived the integrable Korteweg-de Vries equation from the
continuum limit of the FPUT lattice. Introducing energy into one normal mode
with wave number $k$ is nothing else but taking the sinusoidal initial
condition in real space. As time evolves, the state breaks into a series of
localized solutions, i.e., solitons that move and interact with the fixed
ends, i.e. boundaries. Upon interacting with the fixed ends, the solitons
bounce back and return to their initial positions, i.e., giving rise to FPUT
recurrences. Another explanation to the inefficient energy transfer among
normal modes was provided by Izrailev and Chirikov [10] who used the concept
of the overlap of nonlinear resonances. They associated equipartition of
energy with dynamical chaos and were able to estimate a threshold that
separates regular from chaotic dynamics.

Another direction in the study of the FPUT lattice is that of heat
conductivity in the presence of disorder. The main interest is in its
interplay with nonlinearity. For harmonic disordered systems, all eigenmodes
of the infinite system, i.e., Anderson modes, are known to be localized and
form a complete basis [11]. As a linear superposition of Anderson modes, an
initially localized wave in the infinite chain will remain localized at any
time. Whether this behavior changes qualitatively by the introduction of
nonlinearity is still an open question (see for example [12, 13, 14] and
references therein). Disorder can be introduced in the form of uniformly
distributed random variation of particle masses [15, 12], linear coupling
constants between nearest neighbours [16], or in the nonlinearity coefficients
[17]. Recently, by viewing FPUT lattices as systems of masses coupled with
nonlinear springs, Nelson et al. [18] incorporated heterogeneity on a one-
dimensional FPUT array to take into account uncertainties (i.e., in the
masses, the spring constants, or the nonlinear coefficients) during the
manufacturing process of such physical systems. They demonstrated numerically
that tolerances degrade the observance of recurrences, often leading to a
complete loss in moderately-sized arrays. Such a variability may therefore
provide a plausible explanation to little experimental evidence on FPUT energy
recurrences.

Here, we consider the problem of heterogeneous FPUT systems studied in [18].
In our work, we perform numerical simulations in great detail to understand
the breakdown of FPUT recurrences in the model. Indeed, we observe recurrence
degradation, where the energy peak of the lowest normal mode is decreasing
subsequently. For percentage of variability smaller than a threshold that we
derive, the energy is then localized in the few lowest normal modes. The
authors in [19, 20] considered non-equipartition of energy among normal modes
and studied time-periodic states that are exponentially localized in the $k$\-
(or $q$\- in [19, 20]) space of normal modes. Such time-periodic states are
referred to as $q$-breathers. Variability in FPUT lattices therefore leads to
$q$-breathers. In our work, by transforming the FPUT system into another
system in the normal-mode space and considering a two normal-mode
approximation, we provide a qualitative explanation for the disappearance of
FPUT recurrences. In this approximation, $q$-breathers are periodic solutions
centered around an equilibrium point (i.e., time-independent solutions) in the
$q$-space of normal modes.

We also perform long-term numerical integrations and compute the maximum
Lyapunov exponent (mLE) [21] and Smaller Alignment Index (SALI) [22, 23] to
show that the trajectories of the heterogeneous system become quickly chaotic,
as the number of particles in the system increases for the same percentage of
variability. In homogeneous FPUT lattices (i.e., in the absence of
variability), when recurrences occur, the system reaches a metastable state
[24, 25], where only few (low $k$) normal modes share the total energy of the
system. However, it has also been shown that a rather weak diffusion takes
place in the highest normal modes of the spectrum [26] that gradually leads to
equipartition of energy [27]. Using the result in [10], this weak diffusion
process implies weak chaos [28]. Our work shows that variability enhances the
chaotic dynamics of the system.

In this work, we also show that for percentages of variability bigger than a
threshold, solutions may blow up in finite time. Using the same two normal-
mode approximation, we have been able to explain the blow up phenomenon. A
bifurcation analysis is further provided that yields a variability threshold
for the blow up of solutions.

The paper is organised as follows: In Sec. 2, we review the original FPUT
lattice with a quadratic nonlinearity (i.e., the FPUT$-\alpha$ system) and
discuss energy recurrence. We introduce the governing equations of motion in
the presence of parameter variability in Sec. 3. The phenomena of recurrence
breakdown and blow up of solutions are reported in the same section. In Sec.
4, a two normal-mode approximation in the normal mode space is derived using
multiple-scale analysis. Our analytical results explain why energy recurrences
breakdown when variability is introduced and provide a qualitative reason why
solutions blow up in finite-time after a variability threshold. In Sec. 5, we
discuss chaos in the FPUT$-\alpha$ system with or without variability and the
mLE and SALI methods that we use to discriminate between ordered and chaotic
trajectories. Finally, we conclude our study and discuss future work in Sec.
6.

##  2 Mathematical model and dynamics of FPUT$-\alpha$ lattices

The Hamiltonian of the FPUT$-\alpha$ system is given by

$H(x,p)=\frac{1}{2}\sum_{j=0}^{N}p_{j}^{2}+\sum_{j=0}^{N}\frac{1}{2}\left(x_{j+1}-x_{j}\right)^{2}+\frac{\alpha}{3}\left(x_{j+1}-x_{j}\right)^{3}=E,$
(1)

where fixed boundary conditions $x_{0}=x_{N+1}=0$ and $p_{0}=0$ are
considered. In this context, $\alpha\geq 0$ is the nonlinear coupling strength
and $E$ the total, fixed, energy of the system. By viewing the FPUT$-\alpha$
lattice as a model of particles coupled with springs, $x_{j}(t)$ represents
the relative displacement of the $j$th-particle from its equilibrium position
at any time $t$ and $p_{j}(t)$ its corresponding conjugate momentum at any
time $t$. The equations of motion that result from Hamiltonian (1) (i.e.,
Hamilton’s equations of motion) are then given by

$\displaystyle\ddot{x}_{j}=$
$\displaystyle(x_{j+1}-x_{j})+\alpha(x_{j+1}-x_{j})^{2}-(x_{j}-x_{j-1})-\alpha(x_{j}-x_{j-1})^{2}.$
(2)

Working in the real space $x_{j}$ and $p_{j}$, one can express Eqs. (2) in the
normal-mode space $Q_{j}$ and $P_{j}$. This can be done by writing the
position $x_{j}(t)$ as a superposition of eigenvectors of the linear equation.
Using the normal mode transformation,

$\displaystyle\mathbf{x}$
$\displaystyle=A\mathbf{Q},\quad\mathbf{p}=A\mathbf{P},$ (3)

where $\mathbf{x}=[x_{1}~{}x_{2}~{}\ldots~{}x_{N}]^{T}$,
$\mathbf{p}=[p_{1}~{}p_{2}~{}\ldots~{}p_{N}]^{T}$,
$\mathbf{Q}=[Q_{1}~{}Q_{2}~{}\ldots~{}Q_{N}]^{T}$,
$\mathbf{P}=[P_{1}~{}P_{2}~{}\ldots~{}P_{N}]^{T}$, and

$A=\sqrt{\frac{2}{N+1}}\begin{bmatrix}\sin\left(\frac{\pi}{N+1}\right)&\sin\left(\frac{2\pi}{N+1}\right)&\dots&\sin\left(\frac{N\pi}{N+1}\right)\\\
\sin\left(\frac{2\pi}{N+1}\right)&\sin\left(\frac{4\pi}{N+1}\right)&\dots&\sin\left(\frac{2N\pi}{N+1}\right)\\\
\vdots&\vdots&\ddots&\vdots\\\
\sin\left(\frac{N\pi}{N+1}\right)&\sin\left(\frac{2N\pi}{N+1}\right)&\dots&\sin\left(\frac{N^{2}\pi}{N+1}\right)\end{bmatrix},$
(4)

the Hamiltonian (1) becomes

$\displaystyle H$
$\displaystyle=\frac{1}{2}\sum_{k=1}^{N}\left(P_{k}^{2}+\omega_{k}^{2}Q_{k}^{2}\right)+\alpha
H_{3}(Q_{1},Q_{2},\ldots,Q_{N}),$

for some nonlinear function $H_{3}$, where

$\omega_{k}=2\sin\left(\frac{k\pi}{2(N+1)}\right).$ (5)

In this framework, $\mathbf{Q}$ represents the amplitude of the normal mode,
while $\mathbf{P}$ its velocity. The energy of normal mode $k$ for $\alpha=0$
can then be defined by

$\displaystyle E_{k}$
$\displaystyle=\frac{1}{2}\left(P_{k}^{2}+\omega_{k}^{2}Q_{k}^{2}\right).$ (6)

Substituting Eq. (3) into Eq. (2), we obtain the equations of motion in
normal-mode coordinates as

$\displaystyle\ddot{\mathbf{Q}}$
$\displaystyle=D\mathbf{Q}+A^{-1}\mathbf{F(Q)},$ (7)

where

$\displaystyle D=\begin{bmatrix}-\omega_{1}^{2}&0&\ldots&0\\\
0&-\omega_{2}^{2}&&0\\\ \vdots&&\ddots&\vdots\\\
0&0&\ldots&-\omega_{N}^{2}\end{bmatrix},\;\;\mathbf{F(Q)}=\begin{bmatrix}f_{1}(\mathbf{Q})\\\
f_{2}(\mathbf{Q})\\\ \vdots\\\ f_{N}(\mathbf{Q})\end{bmatrix}$

and $A^{-1}$ is the inverse matrix of $A$, given by Eq. (4).

In their seminal work, Fermi, Pasta, Ulam and Tsingou [3] excited the lowest
possible normal mode, i.e., the mode with $k=1$. The initial conditions of
Eqs. (2) are then

$\displaystyle p_{j}$ $\displaystyle=0,\quad x_{j}=\sin\left(\frac{\pi
i}{N+1}\right),\;j=1,2,\ldots,N,$ (8)

which are equivalent to solving Eq. (7) with $Q_{1}=\sqrt{(N+1)/2}$, $Q_{k}=0$
for $k=2,3,\dots,N,$ and $\dot{Q}_{k}=0$ for $k=1,2,\dots,N$.

In Fig. 1, we plot the dynamics of $x_{j}$ of system (1) for $N=64$, for the
initial condition (8), where $E=0.03795$. Panels (a) and (b) show the dynamics
of $x_{j}(t)$ in real space (infact, what is shown is the oscillation
envelope) and the normal mode energy of the first four normal modes of the
FPUT$-\alpha$ system (1), respectively.

(a)

(b)

Figure 1: FPUT recurrences of the Hamiltonian system (1). (a) Dynamics of
$x_{j}(t)$ using the initial condition in (8), where $E=0.03795$. Panel (a)
shows the top view of the oscillation envelope of $x_{j}(t)$ in time. (b)
Energy of the first four normal modes in the dynamics shown in panel (a). Note
in panel (b) how almost all of the energy returns to the first normal mode at
around $t=6\times 10^{4}$, i.e., the appearance of an FPUT recurrence. Here we
have used $N=64$ in the computations in both panels. The range of values in
the vertical axis in panel (a) is between 1 and $N=64$.

In their seminal paper [3], Fermi, Pasta, Ulam and Tsingou expected that the
energy $E$, which was initially used to excite the lowest normal mode only
(i.e., $k=1$), would slowly drift to the other normal modes until the system
reaches thermalization, as predicted by Statistical Mechanics. Surprisingly,
the numerical experiment showed that that was not the case and that after
several periods of the evolution of the mode, almost all energy in the system
returned to the first normal mode that was excited initially. The authors
witnessed the so-called FPUT recurrences. An example of such recurrences is
given in Fig. 1 for $N=64$ and $E=0.03795$.

##  3 Disordered FPUT lattices

The authors in [18] proposed various disordered FPUT$-\alpha$ systems that
include tolerances $v$ into each particle as a result of variability in a
manufacturing process. In their study, they claimed to incorporate tolerances
into the system in different ways based on manufacturing constraints. They
first introduced the following system with heterogeneity

$\displaystyle\ddot{x}_{j}=$
$\displaystyle(v_{j+1}x_{j+1}-v_{j}x_{j})+\alpha(v_{j+1}x_{j+1}-v_{j}x_{j})^{2}-(v_{j}x_{j}-v_{j-1}x_{j-1})-\alpha(v_{j}x_{j}-v_{j-1}x_{j-1})^{2},$
(9)

where, again $\alpha\geq 0$ is the nonlinear coupling strength. It can be
shown that this system admits the Hamiltonian function

$H(x,p)=\frac{1}{2}\sum_{j=0}^{N}\frac{p_{j}^{2}}{v_{j}}+\sum_{j=0}^{N}\frac{1}{2}\left(v_{j+1}x_{j+1}-v_{j}x_{j}\right)^{2}+\frac{\alpha}{3}\left(v_{j+1}x_{j+1}-v_{j}x_{j}\right)^{3}=E.$
(10)

The disorder is inserted in a symmetric way between the linear and nonlinear
coupling. Particularly, the variabilities $v_{j}$ were generated randomly from
a Gaussian distribution, that is for a tolerance $\tau\%$, the values of
$v_{j}$ were drawn from a Gaussian distribution with mean 1 and standard
deviation $\sigma=1/3\times 0.01\tau$. Therefore, the values of $v_{j}$ would
lie in the interval $[1-0.01\tau,1+0.01\tau]$.

The authors in [18] considered also the case where the variabilities $v_{j}$
are present only in the nonlinear coupling terms, resulting in the following
system of second-order ordinary differential equations

$\displaystyle\ddot{x}_{j}=$
$\displaystyle(x_{j+1}-x_{j})+\alpha(v_{j+1}x_{j+1}-v_{j}x_{j})^{2}-(x_{j}-x_{j-1})-\alpha(v_{j}x_{j}-v_{j-1}x_{j-1})^{2}.$
(11)

In this case, the system is no longer Hamiltonian. Then, they showed
numerically that incorporating variability in the nonlinear coupling terms
only has, for a fixed amount of variability, a comparable effect to
incorporating it in only the linear coupling terms. Although in both setups
recurrences such as those in Fig. 1 disappear for large enough tolerance and
the energy localizes in the first few normal modes, more energy is transferred
to the lower modes in the latter case than in the former.

In this work, we consider the effect of disorder in the second scenario of Eq.
(11), which is a toy-dyamical system that does not necessarily relate to a
real physical system. We have decided to study it as we show in Sec. 4, we can
derive a mathematical theory to understand the effect of variability in the
localization of energy and its dynamics.

Throughout the paper, we consider disorder that is generated using the same
setup as in [18]. We show in Fig. 2 the dynamics of $x_{j}$ and normal mode
energies $E_{k}$ of the first four modes of the system of equations (11)
(i.e., for $k=1,\ldots,4$) for $N=64$ particles and two different percentages
of tolerance, i.e., for $\tau=5\%$ and $10\%$.

(a)

(b)

(c)

(d)

Figure 2: Energy recurrences in time for the system in Eq. (11) for $N=64$,
similarly to Fig. 1 for the Hamiltonian system (1). Panels (a) and (b) are for
$\tau=5\%$ tolerance and panels (c) and (d) for $\tau=10\%$. Note that the
ranges in the vertical axes in panels (a) and (c) are from 1 to $N=64$.

Comparing panels (a) in Figs. 1 and 2, one can see that the variability
reduces the effectiveness of recurrence, where a subsequent peak of the mode
energy $E_{1}$ is lower than the preceding ones. In [18], it was reported that
for larger variability, the energy transfer from the lowest to the higher ones
becomes ineffective, which creates a non-recurrent state, shown in panel (d)
in Fig. 2. This state is localized in the normal mode space, i.e., it is a
$q$-breather [19, 20]. In other words, disorder promotes the occurrence of
$q$-breathers.

In Sec. 4, applying a two normal-mode approximation to Eqs. (2) and using
multiple-scale expansions, we show that there is a threshold for the
percentage of variability $\tau_{c}\approx 10.0749\%$, after which the initial
condition (8) may lead to finite-time blow up of the solutions. This is
illustrated in Fig. 3, where the effect of the blow up is clearly seen in the
abrupt increase of the normal-mode energies $E_{k}$ (see Eq. (6)) of the first
four normal modes (i.e., for $k=1,\ldots,4$) for $\tau=20\%>\tau_{c}$.

Figure 3: A similar simulation as in panels (b) and (d) in Fig. 2, for the
system in Eq. (11) and $N=64$, but for the increased tolerance $\tau=20\%$,
where a finite-time blow up of the solution manifests as the abrupt increase
of the energy of the first four normal modes around $t=2000$.

In the following, we show why the transfer of energy between modes reduces
with the increase of the percentage of variability. This leads to a localized
state in the energy-mode space and to finite-time blow up of solutions if
variability is greater than $\tau_{c}$. When energy localization occurs, the
plots of the normal mode energies suggest that most of the mode coordinates
are vanishing in time. Therefore, we prefer to work in the normal-mode
coordinate system than in the real (physical) space. In this framework, the
equations of motion (11) can be written in the normal-mode coordinates space
in a similar manner as in Eq. (7), namely in the form

$\displaystyle\ddot{\mathbf{Q}}$
$\displaystyle=D\mathbf{Q}+A^{-1}\mathbf{\widehat{F}(Q)},$ (12)

for some nonlinear, vector-function $\mathbf{\widehat{F}(Q)}$ that depends on
$\tau$, which is different to $\mathbf{F(Q)}$ in Eq. (7) in the absence of
variability. Our main assumption is that we can approximate system (12) by
considering only the first few modes. To illustrate numerically that this
assumption is reasonable, we present in Fig. 4 the normal-mode energy for the
set of equations of motion (12) for 2, 4 and 8 normal modes and different
percentages of variability. The parameter values in the set of equations of
motion (12) are calculated numerically for $N=64$ and the same percentage of
variability as in Figs. 1 and 2, where all remaining modes are set to 0 at all
times.

Particularly, looking at Fig. 4, we see that using $2$ and $4$ modes gives
dynamics that are quantitatively different from those in Figs. 1 and 2, with
respect to the recurrence period. Nevertheless, even with only $2$ modes, we
can still observe energy recurrence and localization for increasing percentage
of variability. Therefore, in the following, we will consider a two normal-
mode system in Eq. (12).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 4: Normal-mode energy in time obtained from integrating Eq. (12) using
only 2 normal modes in panels (a), (d), (g), 4 modes in panels (b), (e), (h)
and 8 modes in panels (c), (f), (i). The tolerance is $0\%$ in panels (a) -
(c), $5\%$ in panels (d) - (f), and $10\%$ in panels (g) - (i). We note that
for illustration purposes, we plot in all panels only the normal mode energy
of the first four modes and that all tolerances are smaller than $\tau_{c}$.
Despite the fact that the last four modes are activated for $0\%$ and $5\%$
tolerance, they are essentially zero for $10\%$ tolerance in panel (i).

##  4 A two normal-mode system and bifurcation analysis

Figure 4 suggests that when energy localization in the first few normal mode
occurs, all higher modes have relatively much smaller energy. This gives us
the idea that we can approximate Eq. (12) by setting $Q_{k}(t)=0$ for
$k=3,4,\dots,N$, and obtain the following two normal-mode system

$\displaystyle\ddot{Q}_{1}$
$\displaystyle=-\omega_{1}^{2}Q_{1}+\epsilon\left(A_{1}Q_{1}^{2}+A_{2}Q_{2}^{2}+A_{3}Q_{1}Q_{2}\right),$
(13a) $\displaystyle\ddot{Q}_{2}$
$\displaystyle=-\omega_{2}^{2}Q_{2}+\epsilon\left(B_{1}Q_{1}^{2}+B_{2}Q_{2}^{2}+B_{3}Q_{1}Q_{2}\right),$
(13b)

where $A_{i}$, $B_{i}\in\mathbb{R}$, $i=1,2,3$ and $\omega_{k}$ is given in
Eq. (5).

###  4.1 Multiple-scale expansions

Since $\omega_{2}=2\omega_{1}+\epsilon$, $|\epsilon|\ll 1$, we take the
following asymptotic series

$\displaystyle Q_{1}$ $\displaystyle=X_{0}(t,T)+\epsilon X_{1}(t,T)+\ldots,$
(14a) $\displaystyle Q_{2}$ $\displaystyle=Y_{0}(t,T)+\epsilon
Y_{1}(t,T)+\ldots,$ (14b)

where $T=\epsilon t$ is a slow-time variable. The leading-order approximations
to Eqs. (14) are given by

$\displaystyle X_{0}$
$\displaystyle=q_{1}(T)e^{i\omega_{1}t}+q_{1}^{*}(T)e^{-i\omega_{1}t},\quad
Y_{0}=q_{2}(T)e^{i\omega_{2}t}+q_{2}^{*}(T)e^{-i\omega_{2}t}.$ (15)

Substituting Eqs. (14), (15) into Eq. (13), expanding the equations in
$\epsilon$ and applying the standard solvability condition to avoid secular
terms appearing (see e.g., [29]), we obtain

$\displaystyle i\frac{dq_{1}(T)}{dT}$
$\displaystyle=q_{1}(T)+\widetilde{A}q_{1}^{*}q_{2},$ (16a)
$\displaystyle\quad i\frac{dq_{2}(T)}{dT}$
$\displaystyle=q_{2}(T)+\widetilde{B}q_{1}^{2},$ (16b)

for $q_{1}$ and $q_{2}$, respectively, where $\tilde{A}=A_{3}/(2\omega_{1})$
and $\tilde{B}=B_{1}/(2\omega_{2})$. In this context, $i$ is the imaginary
unit of the complex numbers. Following Eqs. (8), the initial conditions of
system (16) are given by

$\displaystyle q_{1}(0)$ $\displaystyle=\frac{Q_{1}(0)}{2},$ (17)
$\displaystyle q_{2}(0)$ $\displaystyle=0.$ (18)

We note that parameters $\tilde{A}$ and $\tilde{B}$ depend on $\tau$.

In Fig. 5, we plot these parameters as a function of $\tau$ for $N=64$
particles and 100 realizations. These realizations have been computed by
fixing $\tau$ and then opting for 100 sets of $N=64$ randomly generated
numbers from the Gaussian distribution with mean 1 and standard deviation
$\sigma=1/3\times 0.01\tau$. Therefore, the $v_{j}$s in the 100 sets lie in
the interval $[1-0.01\tau,1+0.01\tau]$. As we can see in panel (a),
$\tilde{A}$ is positive for all $\tau$, whereas $\tilde{B}$ changes sign at
around $\tau=10\%$. Particularly, $\tilde{B}$ starts positive for small $\tau$
values before it becomes negative at around $\tau=10\%$. By using polynomial
regression, we have been able to fit the mean of the 100 realisations in panel
(b) by the function $\tilde{B}\approx-0.00893\tau^{2}-0.000084\tau+0.90728$,
with a sum of square errors (SSE) of $3.46\times 10^{-19}$. This allowed us to
estimate with good accuracy the threshold for the percentage of variability
where $\tilde{B}$ changes sign and found to be given by $\tau_{c}\approx
10.0749\%$ as $\tilde{B}(\tau_{c})=0$. In Sec. 4.2, we show that when
$\tilde{B}<0$, that is for $\tau>\tau_{c}$, trajectories of Eqs. (2) may blow
up in finite time.

(a)

(b)

Figure 5: Plot of $\tilde{A}$ (in panel (a)) and $\tilde{B}$ (in panel (b)) as
a function of the tolerance obtained numerically for $N=64$. The dash-dotted
curve is the mean value over $100$ realisations of the same percentage of
variability (see the discussion in the text), while the lengths of the shaded
regions are two standard deviations. Using a polynomial regression, the mean
is found to be given approximately by $\tilde{A}\approx
0.01739\tau^{2}-0.00029\tau+3.62805$ and
$\tilde{B}\approx-0.00893\tau^{2}-0.000084\tau+0.90728$, where the sums of
square errors are $2.14\times 10^{-15}$ and $3.46\times 10^{-19}$ in panels
(a) and (b), respectively. Note the horizontal black dashed line at
$\tilde{B}=0$ from which $\tau_{c}$ is derived (see text for details).

A comparison of the dynamics of the normal modes $Q_{1}$ and $Q_{2}$ of Eq.
(13) and those of the slow-time variables $q_{1}$ and $q_{2}$ of Eqs. (16) is
shown in Fig. 6, where one can see that $q_{j}$ is an envelope of $Q_{j}$ for
$j=1,2$.

(a)

(b)

Figure 6: Time evolution of the normal mode variables $Q_{1}$ (blue curve) and
$Q_{2}$ (red curve) with their envelopes $q_{1}$ and $q_{2}$ (black curves)
from Eqs. (16) for $\tau=0\%$ in panel (a) and $\tau=10\%$ in panel (b). Note
that in both panels $\tau<\tau_{c}$, so trajectories do not blow up.

Next we explain the cause of localization with the increase of the percentage
of variability $\tau$. Note that from Eqs. (16), there can be transfer of
energy from $q_{1}(t)$ to $q_{2}(t)$ through the nonlinear coupling
coefficient $\tilde{B}$. Panel (b) in Fig. 5 shows that $\tilde{B}$ decreases
from positive values with the increase of $\tau$ until $\tau=\tau_{c}$, after
which it becomes negative. When $\tilde{B}$ vanishes at $\tau=\tau_{c}$, there
is no transfer of energy and hence localization. In the following, we will
also show that when $\tilde{B}<0$, i.e., for $\tau>\tau_{c}$, there might be
unbounded trajectories that blow up in finite time.

###  4.2 Equilibrium solutions

We start by analyzing the standing wave solutions of the envelope equations
(16). To do so, it is convenient to write $q_{1}$ and $q_{2}$ in polar form
$q_{1}=r_{1}e^{i\phi_{1}}$ and $q_{2}=r_{2}e^{2i\phi_{2}}$, where
$r_{1}=|q_{1}|$, $r_{2}=|q_{2}|$. Then, we define the new variables

$\displaystyle P$ $\displaystyle=r_{1}^{2}+r_{2}^{2},$ (19a)
$\displaystyle\Delta$ $\displaystyle=r_{1}^{2}-r_{2}^{2},$ (19b)
$\displaystyle\theta$ $\displaystyle=\phi_{2}-\phi_{1}.$ (19c)

These variables satisfy the set of equations (see [30] for a similar
derivation)

$\displaystyle\dot{P}$
$\displaystyle=\frac{\widetilde{A}-\widetilde{B}}{\widetilde{A}+\widetilde{B}}\dot{\Delta},$
(20a) $\displaystyle\dot{\Delta}$
$\displaystyle=\frac{\sqrt{2\left(P-\Delta\right)}\sin\left(2\,\theta\right)\left(P+\Delta\right)\left(\widetilde{A}+\widetilde{B}\right)}{2},$
(20b) $\displaystyle\dot{\theta}$
$\displaystyle=-{\frac{2\,\widetilde{A}\cos\left(2\,\theta\right)(\Delta-P)+\widetilde{B}\cos\left(2\,\theta\right)(\Delta+P)-\sqrt{2\left(P-\Delta\right)}}{2\sqrt{2\left(P-\Delta\right)}}}.$
(20c)

From Eq. (20a), the constant of motion $C$ follows

$\displaystyle C$
$\displaystyle=P-\frac{\widetilde{A}-\widetilde{B}}{\widetilde{A}+\widetilde{B}}\Delta.$

System (16) is transformed into Eqs. (20) by using Eqs. (19), where we assume
that $r_{1}$ and $r_{2}$ are non negative real numbers. Equation (20) requires
$P-\Delta>0$ in order to have real-valued solutions, whereas Eq. (19) requires
$P-\Delta>0$ and $P+\Delta\geq 0$, otherwise $r_{1}$ and $r_{2}$ will be
complex numbers. We call the region which satisfies these two inequalities the
well-defined region and denote it by the shaded area in Fig. 7. $\Delta_{2}$
is outside the shaded region in the area below the red curve and above
$\widetilde{B}=0$. This implies that $\Delta_{2}$ is the equilibrium of system
(20) only, but not of system (16).

Figure 7: Bifurcation diagram of the equilibrium points $\Delta_{1}$ and
$\Delta_{2}$ and the regions where the dynamics of system (20) is well-defined
(see text for more details).

The latter result implies that the dynamics of Eq. (16) can be described by
the remaining equations (20b) and (20c), i.e., in terms of $\Delta$ and
$\theta$ only. As discussed before, Eqs. (20b), (20c) are valid only when
$P-\Delta>0$. However, Eqs. (19a), (19b) imply that $P+\Delta\geq 0$. These
two inequalities determine the region where Eqs. (19a), (19b) are defined in
the $(\Delta,\theta)$-plane. As this region depends on $\widetilde{A}$ and
$\widetilde{B}$, we consider the following cases:

  * 1.

If $\frac{\widetilde{A}-\widetilde{B}}{\widetilde{A}+\widetilde{B}}\geq 1$,
then
$\Delta>\text{max}\left\\{\Delta_{\text{crit}}^{1},\Delta_{\text{crit}}^{2}\right\\}$,
where

$\displaystyle\Delta_{\text{crit}}^{1}=\frac{C\left(\widetilde{A}+\widetilde{B}\right)}{2\widetilde{B}}\mbox{
and
}\Delta_{\text{crit}}^{2}=-\frac{C\left(\widetilde{A}+\widetilde{B}\right)}{2\widetilde{A}}.$

  * 2.

If $-1\leq\frac{\widetilde{A}-\widetilde{B}}{\widetilde{A}+\widetilde{B}}<1$,
then $\Delta_{\text{crit}}^{2}\leq\Delta<\Delta_{\text{crit}}^{1}$.

  * 3.

If $\frac{\widetilde{A}-\widetilde{B}}{\widetilde{A}+\widetilde{B}}<-1$, then
$\Delta<\text{min}\left\\{\Delta_{\text{crit}}^{1},\Delta_{\text{crit}}^{2}\right\\}$.

The regions, in which the reduced system (20) is well-defined, are plotted in
the ($\tilde{A}$,$\tilde{B}$)-plane in Fig. 7.

To study the reduced system of Eqs. (20b), (20c), we restrict the phase
difference $\theta$ in the interval $0\leq\theta<\pi$ and obtain two
equilibrium points, namely $(\theta_{j},\Delta_{j})$, $j=1,2$, where

$\theta_{1}=0\mbox{ or }\pi/2\mbox{ and }\theta_{2}=\pi/2.$ (21)

Particularly, there are two cases with respect to $\tilde{B}$. The first one
is when $\widetilde{B}>0$, in which case $\theta_{1}=0$ or $\theta_{2}=\pi/2$,
and the second when $\widetilde{B}<0$, in which case
$\theta_{1}=\theta_{2}=\pi/2$. Then,

$\displaystyle\Delta_{j}$
$\displaystyle={\frac{\left(6\,{\widetilde{A}}^{2}C-3\,\widetilde{A}\widetilde{B}C-(-1)^{j}\sqrt{1+6\,\widetilde{A}\left(\widetilde{A}+\widetilde{B}\right)C}-1\right)\left(\widetilde{A}+\widetilde{B}\right)}{18{\widetilde{A}}^{2}\widetilde{B}}}.$
(22)

The stability of the equilibrium points is determined by the eigenvalues of
the Jacobian matrix of Eqs. (20b), (20c), evaluated at the equilibrium points,
i.e., by

$\displaystyle\lambda^{(j)}_{1,2}$
$\displaystyle=\pm\frac{\sqrt{-3-18{\widetilde{A}}^{2}C-18\widetilde{A}\widetilde{B}C+6(-1)^{j}\sqrt{1+6\,\widetilde{A}\left(\widetilde{A}+\widetilde{B}\right)C}}}{3}.$
(23)

From Eq. (22) it follows that the equilibrium points exist when

$\displaystyle 1+6\,\widetilde{A}\left(\widetilde{A}+\widetilde{B}\right)C\geq
0.$ (24)

For the initial conditions (17), (18), Eq. (24) becomes

$\displaystyle 1+12\widetilde{A}\widetilde{B}r_{1}^{2}\geq 0.$

Therefore, the threshold for the existence of the equilibrium is given by

$\displaystyle 1+12\widetilde{A}\widetilde{B}r_{1}^{2}=0,$

which is the blue curve in Fig. 7. The dashed and solid lines represent the
curve below and above the line $\widetilde{A}+\widetilde{B}=0$, respectively.

Comparing Eqs. (22) and (23), we conclude that when the equilibrium points
exist, they are either a centre or a saddle node. Particularly, for the
initial conditions in Eqs. (17), (18), the thresholds for the eigenvalues that
discriminate between a centre and a saddle node are

$\displaystyle 1+12\widetilde{A}\widetilde{B}r_{1}^{2}$ $\displaystyle=0,$
(25) $\displaystyle 1-4\widetilde{A}\widetilde{B}r_{1}^{2}$ $\displaystyle=0.$
(26)

In Fig. (7), we plot $\tilde{A}+\tilde{B}=0$, Eqs. (25) and (26) as the black
dashed, blue and red curves, respectively.

System (20) with parameter values above the red curve in Fig. (7) is bounded,
with $\Delta_{\text{crit}}^{1}$ and $\Delta_{\text{crit}}^{2}$ being the upper
and lower bounds, respectively. The two equilibria given in Eq. (22) are both
centres, and are therefore stable. When the parameter values lie on the red
curve, $\Delta_{2}=\Delta_{\text{crit}}^{2}$. Furthermore, if the parameter
values are below the red curve and $\widetilde{B}>0$, system (20) is still
bounded, but it only shares one equilibrium point $\Delta_{1}$ with system
(16), whereas $\Delta_{2}$ does not belong to the well-defined region.
Equation (20), on the other hand, is unbounded when $\widetilde{B}<0$. In this
case, it either extends to $\Delta\to\infty$ or $-\infty$ and depending on the
value of $\frac{\widetilde{A}-\widetilde{B}}{\widetilde{A}+\widetilde{B}}$,
$\Delta_{1}$ can be a centre and $\Delta_{2}$ a saddle node in this region.
Additionally, the system has only one equilibrium on the blue curve.

The location of the equilibrium points $(\theta_{j},\Delta_{j})$ in Eqs. (21),
(22) and their nature are shown in Fig. 7. We also plot the values of
$\Delta_{j}$ in Fig. 8. To better visualise $\Delta_{2}$ as it approaches
infinity when $\widetilde{A}$ or $\widetilde{B}$ approaches zero, we plot in
Fig. 8 (b) $\tanh(\Delta_{2}/100)$ instead of $\Delta_{2}$.

(a)

(b)

Figure 8: Plot of (a) $\Delta_{1}$ and (b) $\tanh(\Delta_{2}/100)$ as a
function of $\widetilde{A}$ and $\widetilde{B}$. In panel (a), the color bar
denotes the values of $\Delta_{1}$ and in panel (b), the values of
$\tanh(\Delta_{2}/100)$. The black dashed, red and blue curves are discussed
in the text and are the same with those in Fig. 7.

In the following, we illustrate the phase portrait of the reduced system of
Eqs. (20b), (20c) for different percentages of variability $\tau$, which
correspond to different values of $\tilde{A}$ and $\tilde{B}$. When there is
no variability (i.e., for $\tau=0\%$), the parameter values are
$\widetilde{A}=3.63$ and $\widetilde{B}=0.91$ and the equilibrium points are
$(\theta_{1},\Delta_{1})=(0,5.09)$ and $(\theta_{2},\Delta_{2})=(\pi/2,4.34)$.
Both are stable and the phase space in this case is shown in Fig. 9(a). As we
can see in panel (b) in Fig. 5, as $\tau$ increases, $\widetilde{B}$ decreases
and becomes negative for $\tau>\tau_{c}$. The parameter values for $\tau=10\%$
variability are $\widetilde{A}=4.97$ and $\widetilde{B}=0.05$ and the
equilibrium points are $(\theta_{1},\Delta_{1})=(0,6.27)$ and
$(\theta_{2},\Delta_{2})=(\pi/2,4.08)$. Similar to the previous case, both
equilibrium points are stable and the phase space is shown in Fig. 9(b). Note
that for the initial conditions (18), we have that

$\displaystyle\lim_{\widetilde{B}\rightarrow 0}\Delta_{\text{crit}}^{1}$
$\displaystyle=r_{1}^{2},\quad\lim_{\widetilde{B}\rightarrow
0}\Delta_{\text{crit}}^{2}=0,$

which shows that $\Delta$ becomes positive as we increase $\tau$. Indeed,
$\Delta>0$ corresponds to energy localization as the magnitude of $q_{1}$
remains larger than $q_{2}$.

(a)

(b)

Figure 9: Phase portraits of the reduced system of Eqs. (20b), (20c) for (a)
$\tau=0\%$ percentage of variability and (b) $\tau=10\%$ percentage of
variability.

As we can see in Fig. 10 for $\tau\approx 10.0833\%>\tau_{c}$, $\widetilde{B}$
is negative ($\widetilde{B}=-0.0015$) and the region in the
$(\Delta,\theta)$-space becomes unbounded (see also Fig. 7). It extends to
either $\Delta\to\infty$ or $-\infty$ and depends on $\widetilde{A}$. In this
case, the two equilibrium points are $(\theta_{1},\Delta_{1})=(\pi/2,9.1274)$,
which is a (stable) center, and $(\theta_{2},\Delta_{2})=(\pi/2,15.4383)$,
which is a (unstable) saddle point. The plot shows that in this case, one may
obtain bounded solutions as well as unbounded ones, depending on the initial
condition. For example, the initial condition of Eqs. (8) (or Eqs. (17), (18))
results in $\theta$ and $\Delta$ values in the unbounded region in Fig. 10,
where the trajectory is shown as the blue curve and starts at the bottom of
the plot.

Figure 10: The same as Fig. 9, where the parameter values are
$\widetilde{A}=5.3932$ and $\widetilde{B}=-0.0015$, which correspond to
$\tau\approx 10.0833\%>\tau_{c}$. The blue curve is the trajectory of the
initial condition in Eqs. (17), (18).

##  5 Chaotic behavior

Energy recurrences arise in the homogeneous FPUT lattice (1) when the system
remains in the quasi-stationary state for an extremely long time, making the
approach to equipartition of energy unobservable. In the quasi-stationary
state, the FPUT lattice can be viewed as the perturbation of the regular,
integrable Toda lattice [26].

Here we study the effect of variability on the chaotic properties of system
(11). Particularly, we consider lattices of $N=4,8,16,32,64$ particles in
systems (1) (homogeneous, no variability) and (11) (with variability) and use
the maximum Lyapunov exponent (mLE) [21] and Smaller Alignment Index (SALI)
[22, 23] to discriminate between regular and chaotic dynamics. We want to see
if energy localization in the first normal mode that we observed in Secs. 3
and 4 for $\tau=10\%<\tau_{c}$ corresponds to chaotic dynamics, by increasing
$\tau$ from 0 to $10\%$.

To compute mLE, we follow the evolution of a trajectory starting at the
initial point

$\mathbf{x}(0)=(q_{1}(0),\ldots,q_{N}(0),p_{1}(0),\ldots,p_{N}(0)),$

that evolves according to Hamilton’s equations of motion

$\dot{\mathbf{x}}=\mathbf{f(x)}=\left[\frac{\partial{H}}{\partial{\mathbf{p}}}\quad-\frac{\partial{H}}{\partial{\mathbf{q}}}\right]^{T},$

and the evolution of a deviation vector

$\mathbf{w}(0)=(\delta q_{1}(0),\ldots,\delta q_{N}(0),\delta
p_{1}(0),\ldots,\delta p_{N}(0)),$

that evolves according to the variational equation

$\displaystyle\dot{\mathbf{w}}$
$\displaystyle=\frac{\partial{\mathbf{f}}}{\partial{\mathbf{x}}}(\mathbf{x}(t))\cdot\mathbf{w}.$
(27)

Then mLE is defined as

$\displaystyle\lambda=\lim_{t\to\infty}\frac{1}{t}\ln\frac{||\mathbf{w}(t)||}{||\mathbf{w}(0)||},$

where $\ln$ is the natural logarithm. If mLE converges to zero following the
law $1/t$, then the trajectory is regular, whereas if it converges to a
positive value in time, then the trajectory is chaotic [31]. Hence it is
convenient to plot mLE in $\log_{10}$$-\log_{10}$ scales as the law $1/t$
becomes then a line with negative slope and serves as a guide to the eye.

To compute SALI, we follow the evolution of the same initial condition and two
deviation vectors $\mathbf{w}_{1}(0)$, $\mathbf{w}_{2}(0)$. Then, SALI is
defined by

$\displaystyle\textrm{SALI}(t)=\min\\{\left\|\hat{\mathbf{w}}_{1}(t)-\hat{\mathbf{w}}_{2}(t)\right\|,\left\|\hat{\mathbf{w}}_{1}(t)+\hat{\mathbf{w}}_{2}(t)\right\|\\},$

where
$\hat{\mathbf{w}}_{i}(t)=\frac{\mathbf{w}_{i}(t)}{\left\|\mathbf{w}_{i}(t)\right\|},\;i=1,2$,
are the two normalized deviation vectors at time $t$. SALI approaches zero
exponentially fast in time (as a function of the largest or 2 largest Lyapunov
exponents) for chaotic trajectories and non-zero, positive, values for regular
trajectories [23].

First, we consider the case without variability, that is the FPUT$-\alpha$
system (1). We integrate the equations of motion (2) and its corresponding
variational equations (following Eq. (27)) by using the tangent-map method
[32] and Yoshida’s fourth order symplectic integrator [33]. We have found that
a time step of $0.01$ keeps the relative energy error below $10^{-9}$. In all
our computations, the final integration time is $t=10^{8}$. Here, we use the
same initial condition in Eq. (8) for all $N$. This initial condition then
results in different energies for different $N$, i.e., $E=0.4775$ for $N=4$,
$E=0.2714$ for $N=8$, $E=0.1447$ for $N=16$, $E=0.0747$ for $N=32$, and
$E=0.0379$ for $N=64$. Our results in Fig. 11 show that all trajectories for
$N=4,8,16,32,64$ are regular up to $t=10^{8}$, corroborated by the tendency of
the mLEs to converge to zero following the $1/t$ law and SALI to tend to fixed
positive values, shown in panels (a) and (b), respectively. These results are
in agreement with the fact that energy recurrences in the homogeneous FPUT
lattice (1) arise when it remains in the quasi-stationary state for extremely
long times, making the approach to equipartition of energy unobservable.

(a)

(b)

Figure 11: Plot of mLE (panel a)) and SALI (panel b)) in time for a range of
$N$ values seen in the insets (denoted by different colors) in the absence of
variability, i.e., of the FPUT system (1). Note that all axes are logarithmic.
The black dashed line in panel (a) is the law $1/t$ of regular trajectories to
guide the eye.

Finally, we look at the case of $\tau=10\%<\tau_{c}$, for which we have
observed almost energy localization in the first normal mode in Sec. 3. Since
in this case we only know the equations of motion (11), we integrated them
using the DOP853 integrator [34], an explicit Runge-Kutta method of order $8$
due to Dormand and Prince, to achieve good numerical accuracy. We compute the
chaotic indicators for $30$ realisations of the same percentage of variability
$\tau=10\%$, while keeping the initial conditions fixed for each number of
particles $N$. For $N=4$ and 8, all trajectories in panels (a)-(d) in Fig. 12
appear to be regular up to final integration time $t=10^{8}$, corroborated by
the tendency of the mLEs to converge to zero following the $1/t$ law and SALI
to tend to fixed positive values. However, for $N=16$, two of the 30
trajectories in panels (e), (f) in Fig. 13 are chaotic as their mLEs converge
to positive values at $t=10^{8}$ and their SALI decrease to zero exponentially
fast. Figure 13 shows that there are more chaotic orbits than those for the
smaller values of $N$ in Fig. 12. We show the percentage of chaotic
trajectories (out of the 30 realisations) as a function of $N$ in Fig. 14,
where the increase from $N=4,8,16$ to $N=32,64$ is apparent. These results
suggest that in the case of almost complete energy localization, variability
promotes chaos in the system as the number of particles increases. However,
further studies are required to determine whether the increase is monotone.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 12: Plot of mLE (panels (a), (c), (e)) and SALI (panels (b), (d), (f))
in time for 30 trajectories (denoted by different colors) and $\tau=10\%$ (see
Eq. (11)). Panels (a), (b) are for $N=4$, panels (c), (d) for $N=8$ and panels
(e), (f) for $N=16$. Note that all axes are logarithmic. The black dashed
lines in panels (a), (c), (e) are the law $1/t$ of regular trajectories to
guide the eye.

(a)

(b)

(c)

(d)

Figure 13: Plot of mLE (panels (a), (c)) and SALI (panels (b), (d)) in time
for 30 trajectories (denoted by different colors) and $\tau=10\%$ (see Eq.
(11)). Panels (a), (b) are for $N=32$ and panels (c), (d) for $N=64$. Note
that all axes are logarithmic. The black dashed lines in panels (a), (c) are
the law $1/t$ of regular trajectories to guide the eye. Figure 14: Percentage
of chaotic trajectories as a function of $N$ for $30$ realisations of
variability with $\tau=10\%$. The black-dash line segments connect the black
points and are there to guide the eye.

##  6 Conclusions and discussion

In this paper, we have considered a disordered FPUT$-\alpha$ system with
variations in its parameters (also called variability) to take into account
inherent manufacturing processes. By using a two normal-mode approximation, we
have been able to explain the mechanism for energy localization and blow up of
solutions for percentage of variability bigger than a threshold, that we have
been able to compute using our theory. Moreover, we have also studied the
effect of variability in the chaotic behavior of the system calculating the
maximum Lyapunov exponent and Smaller Alignment Index for a number of
realizations for the same variability percentage that corresponds to energy-
localization. We have found that, when there is almost energy localization, it
is more frequent for the trajectories to be chaotic with the increase of the
number of particles $N$ for the same percentage of variability, smaller than
the threshold.

Finally, while it has been shown previously that variability leads to energy-
recurrence breakdown and energy localization, we have also shown here that by
increasing the percentage of variability beyond a threshold that we determined
using our theory, the solutions of the system may blow up in finite-time. This
is because we have started with the equations of motion without a Hamiltonian
that would allow us to keep the energy of the system constant [18]. The case
of the Hamiltonian model with heterogeneity, cf. Eq. (9), will be studied in a
future publication.

## Credit authorship contribution statement

Zulkarnain: Investigation, Writing – Original Draft. H. Susanto:
Conceptualization, Supervision, Methodology, Writing – review & editing. C.
Antonopoulos: Conceptualization, Supervision, Methodology, Writing – review &
editing.

## Declaration of Competing Interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Acknowledgement

Z is supported by the Ministry of Education, Culture, Research, and Technology
of Indonesia through a PhD scholarship (BPPLN). HS is supported by Khalifa
University through a Faculty Start-Up Grant (No. 8474000351/FSU-2021-011) and
a Competitive Internal Research Awards Grant (No. 8474000413/CIRA-2021-065).
The authors acknowledge the use of the High Performance Computing Facility
(Ceres) and its associated support services at the University of Essex in the
completion of this work. The authors are also grateful to the referees for
their comments and feedback that improved the manuscript.

## References

  * [1] Rudolf Peierls.  Zur kinetischen theorie der wärmeleitung in kristallen.  Annalen der Physik, 395(8):1055–1101, 1929.
  * [2] Peter Debye.  Zur theorie der spezifischen wärmen.  Annalen der Physik, 344(14):789–839, 1912.
  * [3] E Fermi, J Pasta, and S Ulam.  Los Alamos report LA-1940.  Fermi, Collected Papers, 2:977–988, 1955.
  * [4] Joseph Ford.  The Fermi-Pasta-Ulam problem: paradox turns discovery.  Physics Reports, 213(5):271–310, 1992.
  * [5] GP Berman and FM Izrailev.  The Fermi-Pasta-Ulam problem: fifty years of progress.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 15(1):015104, 2005.
  * [6] Thierry Dauxois, Michel Peyrard, and Stefano Ruffo.  The Fermi-Pasta-Ulam “numerical experiment”: history and pedagogical perspectives.  European Journal of Physics, 26(5):S3, 2005.
  * [7] Giovanni Gallavotti.  The Fermi-Pasta-Ulam problem: a status report, volume 728.  Springer, 2007.
  * [8] Mason A Porter, Norman J Zabusky, Bambi Hu, and David K Campbell.  Fermi, Pasta, Ulam and the birth of experimental mathematics: a numerical experiment that Enrico Fermi, John Pasta, and Stanislaw Ulam reported 54 years ago continues to inspire discovery.  American Scientist, 97(3):214–221, 2009.
  * [9] Norman J Zabusky and Martin D Kruskal.  Interaction of “solitons” in a collisionless plasma and the recurrence of initial states.  Physical Review Letters, 15(6):240, 1965.
  * [10] F. M. Izrailev and B. V. Chirikov.  Statistical properties of a nonlinear string.  In Soviet Physics Doklady, volume 11, page 30, 1966.
  * [11] Philip W Anderson.  Absence of diffusion in certain random lattices.  Physical review, 109(5):1492, 1958.
  * [12] Baowen Li, Hong Zhao, and Bambi Hu.  Can disorder induce a finite thermal conductivity in 1D lattices?  Physical Review Letters, 86(1):63, 2001.
  * [13] Abhishek Dhar and Keiji Saito.  Heat conduction in the disordered Fermi-Pasta-Ulam chain.  Physical Review E, 78(6):061136, 2008.
  * [14] Jinyuan Zhu, Yue Liu, and Dahai He.  Effects of interplay between disorder and anharmonicity on heat conduction.  Physical Review E, 103(6):062121, 2021.
  * [15] Daniel N Payton III, Marvin Rich, and William M Visscher.  Lattice thermal conductivity in disordered harmonic and anharmonic crystal models.  Physical Review, 160(3):706, 1967.
  * [16] S Lepri, R Schilling, and S Aubry.  Asymptotic energy profile of a wave packet in disordered chains.  Physical Review E, 82(5):056602, 2010.
  * [17] Arkady Pikovsky.  First and second sound in disordered strongly nonlinear lattices: numerical study.  Journal of Statistical Mechanics: Theory and Experiment, 2015(8):P08007, 2015.
  * [18] Heather Nelson, Mason A. Porter, and Bhaskar Choubey.  Variability in Fermi-Pasta-Ulam-Tsingou arrays can prevent recurrences.  Phys. Rev. E, 98:062210, Dec 2018.
  * [19] S Flach, MV Ivanchenko, and OI Kanakov.  $q$-breathers and the Fermi-Pasta-Ulam problem.  Physical Review Letters, 95(6):064102, 2005.
  * [20] S Flach, MV Ivanchenko, and OI Kanakov.  q-breathers in Fermi-Pasta-Ulam chains: Existence, localization, and stability.  Physical Review E, 73(3):036618, 2006.
  * [21] Giancarlo Benettin, Luigi Galgani, Antonio Giorgilli, and Jean-Marie Strelcyn.  Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them. Part 1: Theory.  Meccanica, 15(1):9–20, 1980.
  * [22] Charalampos Skokos, Chris Antonopoulos, Tassos C Bountis, and Michael N Vrahatis.  How does the Smaller Alignment Index (SALI) distinguish order from chaos?  Progress of Theoretical Physics Supplement, 150:439–443, 2003.
  * [23] Ch Skokos, Ch Antonopoulos, TC Bountis, and MN Vrahatis.  Detecting order and chaos in Hamiltonian systems by the SALI method.  Journal of Physics A: Mathematical and General, 37(24):6269, 2004\.
  * [24] A Ponno and D Bambusi.  Korteweg-de Vries equation and energy sharing in Fermi-Pasta-Ulam.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 15(1):015107, 2005.
  * [25] Dario Bambusi and Antonio Ponno.  On metastability in FPU.  Communications in mathematical physics, 264(2):539–561, 2006.
  * [26] Antonio Ponno, Helen Christodoulidi, Ch Skokos, and Sergej Flach.  The two-stage dynamics in the Fermi-Pasta-Ulam problem: From regular to diffusive behavior.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 21(4):043127, 2011.
  * [27] Luisa Berchialla, Antonio Giorgilli, and Simone Paleari.  Exponentially long times to equipartition in the thermodynamic limit.  Physics Letters A, 321(3):167–172, 2004.
  * [28] Chris G Antonopoulos and Helen Christodoulidi.  Weak chaos detection in the Fermi-Pasta-Ulam-$\alpha$ system using $q$-Gaussian statistics.  International Journal of Bifurcation and Chaos, 21(08):2285–2296, 2011.
  * [29] M Syafwan, H Susanto, and SM Cox.  Discrete solitons in electromechanical resonators.  Physical Review E, 81(2):026207, 2010.
  * [30] J Pickton and H Susanto.  Integrability of PT-symmetric dimers.  Physical Review A, 88(6):063840, 2013.
  * [31] Ch Skokos.  The Lyapunov characteristic exponents and their computation.  In Dynamics of Small Solar System Bodies and Exoplanets, pages 63–135. Springer, 2010.
  * [32] Ch Skokos and Enrico Gerlach.  Numerical integration of variational equations.  Physical Review E, 82(3):036704, 2010.
  * [33] Haruo Yoshida.  Construction of higher order symplectic integrators.  Physics letters A, 150(5-7):262–268, 1990.
  * [34] Nonstiff Differential Equations, DOP853 integrator.  http://www.unige.ch/~hairer/software.html, 2022.