# Energy Efficient Communications in RIS-assisted UAV Networks Based on
Genetic Algorithm

Mohammad Javad-Kalbasi  Department of Electrical Engineering
University of Toronto
Toronto, Canada
<EMAIL_ADDRESS>Mohammed S. Al-Abiad  Department of
Electrical Engineering
University of Toronto
Toronto, Canada
<EMAIL_ADDRESS>Shahrokh Valaee  Department of Electrical
Engineering
University of Toronto
Toronto, Canada
<EMAIL_ADDRESS>
###### Abstract

This paper proposes a solution for energy-efficient communication in
reconfigurable intelligent surface (RIS)-assisted unmanned aerial vehicle
(UAV) networks. The limited battery life of UAVs is a major concern for their
sustainable operation, and RIS has emerged as a promising solution to reducing
the energy consumption of communication systems. The paper formulates the
problem of maximizing the energy efficiency of the network as a mixed integer
non-linear program, in which UAV placement, UAV beamforming, On-Off strategy
of RIS elements, and phase shift of RIS elements are optimized. The proposed
solution utilizes the block coordinate descent approach and a combination of
continuous and binary genetic algorithms. Moreover, for optimizing the UAV
placement, Adam optimizer is used. The simulation results show that the
proposed solution outperforms the existing literature. Specifically, we
compared the proposed method with the successive convex approximation (SCA)
approach for optimizing the phase shift of RIS elements.

###### Index Terms:

RIS-assisted UAV network, energy-efficient communications, mixed-integer non-
linear program, block coordinate descent, genetic algorithm, Adam optimizer.

##  I Introduction

111This work was supported in part by funding from the Innovation for Defence
Excellence and Security (IDEaS) program from the Department of National
Defence (DND).

In recent years, unmanned aerial vehicles (UAVs) have gained immense
popularity due to the high likelihood of establishing line-of-sight
connections with ground nodes, rapid deployment, and adjustable mobility [1],
[2]. With such attributes, UAVs are being employed for various applications
such as surveillance, rescue missions, delivery, and communication, among
others. The development of UAV communication networks has been identified as a
promising solution to address the limited coverage area of existing
terrestrial communication networks, especially in remote areas.

However, the limited battery life of UAVs is a major concern for their
sustainable operation. The conventional UAVs are equipped with communication
systems that consume a significant amount of energy, resulting in limited
flight time. To overcome this challenge, the concept of reconfigurable
intelligent surface (RIS) has emerged as a promising solution to reduce the
energy consumption of communication systems. RIS has the potential to enhance
the signal-to-interference-plus-noise ratio (SINR) of wireless channels, and
thereby reduce the power consumption of the communication system, enabling
UAVs to fly for a longer duration. Moreover, RIS can provide communication
links to the blocked ground users (GUs) that do not maintain direct access to
the UAVs. Therefore, by integrating RISs and UAVs, RIS-assisted UAV networks
reduce power consumption of UAVs while improving coverage and connectivity
significantly [1].

In the current literature, communication-efficient problems were explored for
RIS-assisted UAV networks, e.g., [3]-[15]. In particular, the authors of [3]
explored an adaptive RIS-assisted aerial-terrestrial downlink communication
system between UAVs and multi-users by optimizing RIS element allocation and
reflective coefficients. In [4], the authors investigated UAV-user
communications with RIS assistance to maximize the worst-case secrecy rate,
considering transmitter power allocation, RIS beamforming, and UAV trajectory.
In [5], the authors proposed a RIS-assisted UAV communication system to
maximize the received signal power at GUs by jointly optimizing passive and
active beamforming and UAV trajectory. Moreover, the work in [6] addressed the
energy consumption problem for both orthogonal multiple access (OMA) and non-
orthogonal multiple access (NOMA) cases by jointly optimizing the UAV’s
trajectory and the passive beamforming of the RIS elements.

Several works also considered the use of RIS-assisted systems [7], [8], [9],
which proposed an aerial RIS (ARIS)-assisted system that satisfies the ultra-
reliable low latency communication (URLLC) constraints. In [10], the authors
investigated total transmit power minimization for heterogeneous networks that
use multiple UAVs and RISs. However, they did not consider energy efficiency
which is defined as the ratio of total sum-rate of GUs to the total power
consumption of system. Energy efficiency has been widely used to evaluate the
performance of conventional communication systems that rely solely on ground
based stations [11], [12], [13]. Reference [14] focused on maximizing energy
efficiency for a single ARIS-assisted downlink communication with a single
user and did not consider the multiple ARIS-assisted scenario. In [15], the
authors considered a multiple ARIS configuration to maximize the average
energy efficiency for downlink communication between the base station and the
GU. The majority of previous works focused on optimizing either the sum-rate
or power consumption objective functions. However, most proposed solutions
employed computationally expensive methods such as successive convex
approximation (SCA) and reinforcement learning. Additionally, there is a
higher likelihood of getting trapped in a local minimum when utilizing the SCA
method.

Motivated by the aforementioned limitations of the existing related studies,
this paper focuses on the energy-efficient communication in RIS-assisted UAV
networks. In this context, a RIS is deployed to enhance wireless connectivity
while reducing UAV movement. We formulate the problem of maximizing the energy
efficiency of the network as a mixed integer non-linear program in which UAV
placement, UAV beamforming, On-Off strategy of RIS elements, and phase shift
of RIS elements are simultaneously optimized. Since solving this problem is
computationally intractable, we use the blocked coordinate descent (BCD)
approach. BCD is an optimization algorithm that solves a large-scale
optimization problem by iteratively optimizing a subset of variables while
holding the other variables fixed. We utilize a continuous genetic algorithm
to optimize the phase shift of RIS elements and the beamforming of UAV,
whereas a binary genetic algorithm is employed to optimize the On-Off strategy
of the RIS elements. Moreover, we optimize the placement of UAV by using the
Adam optimizer [16]. Our simulation results show that our proposed method
offers improved energy efficiency as compared to the existing literature.

##  II System Model and Problem Formulation

###  II-A System Model

We consider a RIS-assisted UAV network with one UAV, one RIS, and multiple
ground users (GUs). All GUs and UAV have single antennas. The UAV can provide
good connectivity for the GUs. To further improve performance, the UAV
utilizes a RIS. Fig. 1 shows an example of the considered system model. The
UAV flies at a fixed altitude $Z_{U}$ for $T$ seconds to provide communication
services for $K$ GUs with $0$ altitude. By using a Cartesian coordinate
system, we can represent the horizontal position of the UAV, RIS, and the
$k$-th GU as $W_{U}=~{}{[x_{U},y_{U}]}^{T}$, $W_{R}={[x_{R},y_{R}]}^{T}$, and
$W_{k}={[x_{k},y_{k}]}^{T}$ for $k=1,2,\ldots,K$. We consider that the RIS and
GUs maintain fixed positions, but the position of the UAV can be changed. The
RIS is deployed at the altitude of $Z_{R}$ equipped with a controller and
$M=M_{c}\times M_{r}$ reflecting units (RUs) to form a uniform array (UA).
Specifically, each row of the UA has $M_{r}$ RUs with an equal distance of
$d_{r}$ meters and each column of the UA consists of $M_{c}$ RUs with an equal
distance of $d_{c}$ meters. These RUs can reflect the received signals in
adjustable phase shifts. The phase-shift matrix of the RIS is modeled as the
diagonal matrix
$\mathbf{\Theta}=diag(e^{j\theta_{1}},e^{j\theta_{2}},\ldots,e^{j\theta_{M}})$,
where $\theta_{m}\in[0,2\pi)$ for $m=1,2,\ldots,M$. Considering the power
consumption of RISs due to controlling the phase shift values of the
reflecting elements [15], it is often not energy efficient to turn on all the
RIS elements. Let us define the binary variable $x_{m}\in\\{0,1\\}$, where
$x_{m}=1$ indicates that the $m$-th RIS element is on, and $0$ otherwise.
Defining $\mathbf{\mathcal{X}}=diag(x_{1},\ldots,x_{M})$, the effective phase-
shift matrix of the RIS can be written as
$\mathbf{\Theta}_{e}=\mathbf{\Theta}\mathbf{\mathcal{X}}$.

Figure 1: A typical RIS-assisted UAV communication system.

###  II-B Channel Model

The channel between the UAV and the $k$-th GU is given by [17]

$\mathbf{h}_{k}^{UG}=\sqrt{\frac{\beta_{0}}{({d_{k}^{UG})}^{\alpha_{k}^{UG}}}}\left(\sqrt{\frac{{\kappa}_{k}^{UG}}{{\kappa}_{k}^{UG}+1}}+\sqrt{\frac{1}{{\kappa}_{k}^{UG}+1}}\mathbf{\hat{h}}_{k}^{UG}\right),$

where $d_{k}^{UG}:=\sqrt{{||W_{k}-W_{U}||}^{2}+Z_{U}^{2}}$ is the distance
between the UAV and the $k$-th GU, $\beta_{0}$ is the path loss at the
reference distance of $1$ meter, $\alpha_{k}^{UG}$ represents the path loss
exponent of the link from the UAV to the $k$-th GU, ${\kappa}_{k}^{UG}$
denotes the corresponding Rician factor, and
$\mathbf{\hat{h}}_{k}^{UG}\sim\mathcal{CN}(\mathbf{0},1)$ is the scattering
component of the $k$-th GU.

The link from the UAV to the RIS is presumed by a line-of-sight (LOS) channel,
and the channel gain can be found by [18]

$\mathbf{h}^{UR}=\sqrt{\frac{\beta_{0}}{({d^{UR})}^{2}}}\mathbf{h}_{LOS}^{UR},$

where $d^{UR}:=\sqrt{{||W_{R}-W_{U}||}^{2}+{\left(Z_{R}-Z_{U}\right)}^{2}}$ is
the distance between the UAV and the RIS, and $\mathbf{h}_{LOS}^{UR}$
represents the array response component which can be denoted by

$\displaystyle\mathbf{h}_{LOS}^{UR}=$ $\displaystyle{[1,e^{-j\frac{2\pi
d_{r}}{\lambda}\phi^{UR}\psi^{UR}},\ldots,e^{-j\frac{2\pi
d_{r}}{\lambda}(M_{r}-1)\phi^{UR}\psi^{UR}}]}^{T}$
$\displaystyle\otimes{[1,e^{-j\frac{2\pi
d_{c}}{\lambda}\varphi^{UR}\psi^{UR}},\ldots,e^{-j\frac{2\pi
d_{c}}{\lambda}(M_{c}-1)\varphi^{UR}\psi^{UR}}]}^{T}$

where $\phi^{UR}=\frac{y_{U}-y_{R}}{||W_{R}-W_{U}||}$,
$\varphi^{UR}=\frac{x_{R}-x_{U}}{||W_{R}-W_{U}||}$, and
$\psi^{UR}=~{}\frac{z_{U}-z_{R}}{d^{UR}}$. We point out that $\phi^{UR}$,
$\varphi^{UR}$, and $\psi^{UR}$ are the cosine of horizontal angle-of-arrival
(AoA), sine of horizontal AoA, and sine of vertical AoA, respectively.

Considering the channel between RIS and each GU, we assume that the RIS-aided
channels are modeled by the Rician fading. Thus, the channel gain is given by
[17]

$\mathbf{h}_{k}^{RG}=\sqrt{\frac{\beta_{0}}{({d_{k}^{RG})}^{\alpha_{k}^{RG}}}}\bigl{(}\sqrt{\frac{{\kappa}_{k}^{RG}}{{\kappa}_{k}^{RG}+1}}\mathbf{h}_{k,LOS}^{RG}\\\
+\sqrt{\frac{1}{{\kappa}_{k}^{RG}+1}}\mathbf{\hat{h}}_{k}^{RG}\bigr{)},$

where $d_{k}^{RG}:=\sqrt{{||W_{k}-W_{R}||}^{2}+Z_{R}^{2}}$ is the distance
between the RIS and the $k$-th GU, $\alpha_{k}^{RG}$ denotes the path loss
exponent of the link between the RIS and the $k$-th GU, ${\kappa}_{k}^{RG}$ is
the corresponding Rician factor,
$\mathbf{\hat{h}}_{k}^{RG}\sim\mathcal{CN}(\mathbf{0},I_{M})$ is the
scattering component, and the LOS component is

$\displaystyle\mathbf{h}_{k,LOS}^{RG}=$ $\displaystyle{[1,e^{-j\frac{2\pi
d_{r}}{\lambda}\phi_{k}^{RG}\psi_{k}^{RG}},\ldots,e^{-j\frac{2\pi
d_{r}}{\lambda}(M_{r}-1)\phi_{k}^{RG}\psi_{k}^{RG}}]}^{T}$
$\displaystyle\otimes{[1,e^{-j\frac{2\pi
d_{c}}{\lambda}\varphi_{k}^{RG}\psi_{k}^{RG}},\ldots,e^{-j\frac{2\pi
d_{c}}{\lambda}(M_{c}-1)\varphi_{k}^{RG}\psi_{k}^{RG}}]}^{T},$

in which $\phi_{k}^{RG}=\frac{y_{k}-y_{R}}{||W_{k}-W_{R}||}$,
$\varphi_{k}^{RG}=\frac{x_{k}-x_{R}}{||W_{k}-W_{R}||}$, and
$\psi_{k}^{RG}=~{}\frac{z_{R}}{d_{k}^{RG}}$. In a similar way,
$\phi_{k}^{RG}$, $\varphi_{k}^{RG}$, and $\psi_{k}^{RG}$ denote cosine of
horizontal angle-of-departure (AoD), sine of horizontal AoD, and sine of
vertical AoD, respectively.

Given the aforementioned channel models, the effective channel gain between
the UAV and the $k$-th GU with the aid of the RIS is given by
$\mathbf{C}_{k}=\mathbf{h}_{k}^{UG}+{(\mathbf{h}_{k}^{RG})}^{H}\mathbf{\Theta}_{e}\mathbf{h}^{UR}$.

###  II-C Problem Formulation

In this paper, the UAV is assumed to share the same frequency band for
providing services to the GUs. The transmitted signal of the UAV is given by
$\mathbf{S}=\sum\limits_{k=1}^{K}\sqrt{p_{k}}s_{k}$, where $p_{k}$ and $s_{k}$
are the transmitted power and signal for the $k$-th GU, respectively. Assuming
that $P_{max}$ is the maximum transmit power of the UAV, we have
$\sum\limits_{k=1}^{K}p_{k}\leq P_{max}$. Then the received signal at the
$k$-th GU can be expressed as

$y_{k}=\mathbf{C}_{k}\sqrt{p_{k}}s_{k}+\underbrace{\mathbf{C}_{k}\sum\limits_{t=1,t\neq
k}^{K}\sqrt{p_{t}}s_{t}}_{\text{interference}}+\underbrace{n_{k}}_{\text{noise}},$

where $n_{k}\sim\mathcal{CN}(0,\sigma^{2})$. Therefore, the received signal-
to-interference-plus-noise ratio (SINR) at the $k$-th GU is

$\gamma_{k}=\frac{{|{\mathbf{C}}_{k}|}^{2}p_{k}}{{|{\mathbf{C}}_{k}|}^{2}\sum\limits_{t=1,t\neq
k}^{K}p_{t}+\sigma^{2}}.$

As a result, the sum-rate of all GUs can be written as

$R_{t}=B\sum\limits_{k=1}^{K}\log_{2}(1+\gamma_{k}),$ (1)

where $B$ is the bandwidth of the channel.

Essentially, the UAV consumes power due to hovering and wireless signal
transmission, while GUs consume power for circuit operations and the RIS
elements consume power for controlling the phases. Consequently, the total
power consumption of our envisioned system is given by

$P_{t}=\underbrace{p_{h}}_{\text{hovering power of
UAV}}+\underbrace{\sum\limits_{k=1}^{K}p_{k}}_{\text{transmit power of
UAV}}\\\ +\underbrace{\sum\limits_{k=1}^{K}p_{k}^{c}}_{\text{circuit power of
GUs}}+\underbrace{\sum\limits_{m=1}^{M}p^{r}x_{m}}_{\text{power consumption of
RIS}},$ (2)

where $p_{k}^{c}$ is the circuit power consumption of the $k$-th GU and
$p^{r}$ is the power consumption of each RU. In (2), $p_{h}$ is the drone’s
hovering power which is expressed by $p_{h}=\sqrt{\frac{{{(mg)}}^{3}}{2\pi
r_{p}^{2}n_{p}\rho}}$, where $m$ is the drone’s weight, $g$ is the
gravitational acceleration of the earth, $r_{p}$ is propellers’ radius, $np$
is the number of propellers, and $\rho$ is the air density [19]. From (1) and
(2), we can define the energy efficiency of the RIS-assisted UAV network as
$\eta=\frac{R_{t}}{P_{t}}$.

Given the considered system model, our objective is to jointly optimize the
phase-shift of RIS elements, RIS On-Off vector, beamforming vector of the UAV,
and its placement to maximize the energy efficiency under the minimum rate
requirement of GUs and the total power constraint of the UAV. Mathematically,
the energy efficiency maximization optimization problem can be formulated as

$\displaystyle\underset{\mathbf{X},\boldsymbol{\theta},\mathbf{P},W_{U}}{\text{maximize}}$
$\displaystyle\eta$ subject to $\displaystyle B\log_{2}(1+\gamma_{k})\geq
R_{k,min},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}(3a)$
$\displaystyle\sum\limits_{k=1}^{K}p_{k}\leq
P_{max},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}(3b)$ $\displaystyle
p_{k}>0,~{}~{}~{}~{}~{}\text{for}~{}1\leq k\leq
K,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}(3c)$ $\displaystyle
x_{m}\in\\{0,1\\},~{}~{}~{}~{}~{}\text{for}~{}1\leq m\leq M$
$\displaystyle\theta_{m}\in[0,2\pi),~{}~{}~{}~{}~{}\text{for}~{}1\leq m\leq M$

where $\mathbf{X}={(x_{1},\ldots,x_{M})}^{T}$,
$\boldsymbol{\theta}={(\theta_{1},\ldots,\theta_{M})}^{T}$,
$\boldsymbol{P}=~{}{(p_{1},\ldots,p_{K})}^{T}$, and $R_{k,min}$ is the minimum
data rate requirement of the $k$-th GU. The minimum rate constraint for each
GU is given in (3a). Moreover, (3b) and (3c) represent the total power
constraint of the UAV.

The derived optimization problem is a mixed-integer nonlinear program (MINLP)
even for the single GU case with $K=1$. It is generally difficult to obtain
the globally optimal solution of a MINLP.

##  III Proposed Solution

In this section, we first use the alternating optimization approach (BCD) to
divide the derived MINLP into three sub-problems, which has been widely
applied to tackle the non-concave problems in RIS-assisted UAV networks [18].
Then each sub-problem is separately optimized, while other variables are
considered to be fixed. Specifically, we employ a continuous genetic algorithm
to jointly optimize the phase shifts of the RIS and beamforming of UAV, and a
binary genetic algorithm to optimize the On-Off strategy of the RIS.
Afterwards, the optimized position of the UAV is obtained using the Adam
optimizer. The block diagram of the overall proposed solution is depicted in
Fig. 2. As it can be observed, the optimization process terminates if the
amount of improvement in energy efficiency is less than a given threshold
$\delta$.

UAV placement (Adam optimizer)Phase shift and beamforming (continuous genetic
algorithm)$\left(W_{U}^{(0)},{\mathbf{X}}^{(0)}\right)$On-Off strategy (binary
genetic
algorithm)$W_{U}^{(0)}$$\Delta\eta<\delta$No$\left(\boldsymbol{\theta}^{(i)},\boldsymbol{P}^{(i)}\right)$$W_{U}^{(i)}$${\mathbf{X}}^{(i)}$yes
Figure 2: Block diagram of the overall proposed solution.

Generate the initial population
$\mathbb{P}^{(0)}={[\mathbf{Y}_{1}^{(0)},\mathbf{Y}_{2}^{(0)},\ldots,\mathbf{Y}_{2L}^{(0)}]}^{T}$.

Let $Q^{(0)}$ be a selection random variable defined on the set
$\\{1,2,\ldots,2L\\}$ with probability mass function
$Pr\\{Q^{(0)}=k\\}=\frac{\eta\left({\mathbf{Y}_{k}^{(0)}}\right)}{\sum\limits_{l=1}^{2L}\eta\left({\mathbf{Y}_{l}^{(0)}}\right)}$.

for $i=1:N_{1}$ do

Select $2L$ individuals from $\mathbb{P}^{(i-1)}$ using the selection random
variable $Q^{(i-1)}$.

Randomly select $L$ pairs of individuals in $\mathbb{P}^{(i-1)}$ for cross-
over. Substitute each pair with the two generated offsprings.

Add a small random value to each element of generated offsprings (mutation).

Find the updated population $\mathbb{P}^{(i)}$ and selection random variable
$Q^{(i)}$.

end for

Algorithm 1 Continuous genetic algorithm for phase shift and beamforming
optimization

0: $\alpha$: Step size

0: $\epsilon$: Small constant

0: $\beta_{1},\beta_{2}\in[0,1)$: Exponential decay rates for the moment
estimates

0: $\eta\left(W_{U}\right)$: Objective function

0: $W_{U}^{(0)}$: Initial position of UAV

Initialize $m_{0}\leftarrow 0$, $v_{0}\leftarrow 0$, and $i\leftarrow 0$

for $i=1:N_{2}$ do

Compute gradient
$g_{i}\leftarrow\nabla_{W_{U}}\eta\left(W_{U}^{(i-1)}\right)$.

Update biased first moment estimate:
$m_{i}\leftarrow\beta_{1}m_{i-1}+(1-\beta_{1})g_{i}$

Update biased second raw moment estimate:
$v_{i}\leftarrow\beta_{2}v_{i-1}+(1-\beta_{2})g_{i}^{2}$

Compute bias-corrected first moment estimate:
$\hat{m}_{i}\leftarrow\frac{m_{i}}{1-\beta_{1}^{i}}$

Compute bias-corrected second raw moment estimate:
$\hat{v}_{i}\leftarrow\frac{v_{i}}{1-\beta_{2}^{i}}$

Update parameters: $W_{U}^{(i)}\leftarrow
W_{U}^{(i-1)}-\alpha\frac{\hat{m}_{i}}{\sqrt{\hat{v}_{i}}+\epsilon}$

end for

Algorithm 2 Adam optimizer for the UAV placement. $g_{i}^{2}$ indicates the
element-wise square $g_{i}\odot g_{i}$. Good default settings are
$\alpha=0.001$, $\beta_{1}=0.9$, $\beta_{2}=0.999$, and $\epsilon=10^{-8}$.
All operations on vectors are element-wise. With $\beta_{1}^{i}$ and
$\beta_{2}^{i}$, we denote $\beta_{1}$ and $\beta_{2}$ to the power of $i$.

###  III-A Phase Shift and Beamforming Optimization

For optimizing the phase shift of the RIS elements and beamforming of the UAV,
we take advantage of the continuous genetic algorithm. The continuous genetic
algorithm is a type of optimization algorithm that is used to search for
optimal solutions in continuous parameter spaces [20]. The algorithm is based
on the principles of natural selection and evolution, and it works by creating
a population of candidate solutions and iteratively refining them over
multiple generations. The main steps of the continuous genetic algorithm
include initialization, selection, crossover, mutation, and termination.

  * •

The first step of the algorithm is initialization, where a population of
potential solutions is created randomly. The size of the population is
typically set based on the complexity of the problem, and each candidate
solution is represented as a vector of continuous values.

  * •

In the selection step, the candidate solutions are evaluated based on a
fitness function, which measures how well they perform on the given problem.
In this paper, energy efficiency ($\eta$) is considered as the fitness
function. The best-performing solutions are selected to become parents for the
next generation.

  * •

In the crossover step, pairs of parents are combined to produce new offspring
solutions. This is done by finding the weighted sum of parent vectors and
creating a new solution vector that inherits traits from both parents. The
crossover process helps to explore different regions of the search space and
can help to create more diverse candidate solutions.

  * •

In the mutation step, random changes are introduced into the offspring
solutions to further explore the search space and prevent the algorithm from
getting stuck in a local optimum. This is done by adding a small random value
to each element of the solution vector.

The algorithm continues to iterate through the selection, crossover, and
mutation steps until a termination criterion is met. The details of the
continuous genetic algorithm for optimizing the phase shift of the RIS
elements and beamforming of the UAV are given in Algorithm 1.

###  III-B RIS On-Off Optimization

Given UAV position $W_{U}$, phase vector $\boldsymbol{\theta}$, and
beamforming vector $\boldsymbol{P}$, problem (3) is a nonlinear integer
optimization problem with respect to the RIS On-Off vector $\mathbf{X}$. Since
the nonlinear integer optimization problem is NP-hard in general, it is
difficult to obtain the globally optimal solution with polynomial complexity.
To tackle this computational intractability, we use the binary genetic
algorithm to solve the RIS On-Off optimization problem. The main steps of
binary genetic algorithm are same as the continuous one, where crossover and
mutation steps are adapted to be applied to the vectors of binary values. In
the crossover step, parts of the parent vectors are swapped to create a new
solution vector, whereas binary bits of offspring solutions are flipped with a
small probability in the mutation step.

###  III-C UAV Placement Optimization

To obtain the optimized position of the UAV, we utilize the Adam optimizer.
The Adam optimizer is a gradient-based optimization algorithm that is widely
used in machine learning because of its superior performance compared to other
gradient descent methods [16]. One of the primary reasons for its superiority
is that it maintains separate learning rates for each parameter and adapts
them based on the historical gradient information. This allows the algorithm
to handle sparse gradients and noisy data more efficiently, leading to faster
convergence and better performance. Additionally, the adaptive learning rates
ensure that the algorithm can navigate complex loss landscapes more
effectively, avoiding getting stuck in local minimum.

The Adam optimizer works by computing a moving average of the gradient and its
square, which allows it to adjust the learning rate adaptively. The algorithm
consists of several steps, including initializing the moving averages,
computing the gradient on a mini-batch of data, updating the moving averages
with the gradient information, and finally updating the parameters based on
the adaptive learning rate. The algorithm also incorporates bias correction
terms to ensure that the moving averages are initialized at zero and that the
updates are unbiased. The detailed steps of the Adam optimizer for optimizing
the UAV placement are given in Algorithm 2.

Figure 3: Energy efficiency versus the number of GUs when $M=60$.

##  IV Numerical Results

This section presents numerical results that validate the efficacy of the
proposed alternating approach for optimizing UAV placement and beamforming,
On-Off strategy of RIS elements, and phase shift of RIS elements. The
following benchmark schemes are used for performance comparison:

  * •

Without RIS: This scheme does not employ the RIS to assist communication and
considers only the direct link between the UAV and the GUs.

  * •

RIS with random Phase: In this scheme, the phase shift of the RIS elements are
chosen randomly within the range of $[0,2\pi)$.

  * •

SCA method: In this scheme, we apply SCA method to optimize the phase shift of
the RIS elements. SCA is an iterative optimization technique used to solve
non-convex optimization problems [10]. It involves transforming a non-convex
optimization problem into a sequence of convex sub-problems, which can be
solved efficiently.

In this paper, we examine a scenario similar to the configuration presented in
[18], where a single RIS is installed on a building facade to improve aerial-
ground communications. The GUs are distributed randomly and uniformly within a
circular area centered at $(200~{}m,25~{}m)$ with a radius of $20$ m. The
details of the considered simulation parameters are given in Table I.

Assuming that there are 60 RIS elements, Fig. 3 illustrates the relationship
between the energy efficiency and the number of GUs. The figure shows that the
energy efficiency increases as the number of GUs ($K$) increases.
Additionally, it is apparent that not using RIS results in consistently low
energy efficiency since there is no indirect path from the UAV to the GUs.
Overall, our proposed method is superior to the SCA method and random phase,
particularly for a large number of GUs since the SCA method is more prone to
trapping in a local minimum.

Fig. 4 illustrates the relationship between the energy efficiency and the
number of the RIS elements ($M$) when $K=4$. Fig. 4 clearly shows that
deploying RIS significantly improves the system performance, whereas
negligible gains are observed without the RIS. Optimizing phase shifts results
in slightly better performance compared to the random phase shift scheme,
emphasizing the importance of optimizing phase shifts. The proposed method
outperforms the SCA method, and this improvement becomes more apparent as the
number of GUs increases.

In summary, integrating RIS with UAV communications can significantly enhance
network coverage, and the energy efficiency of the RIS-assisted schemes
improves as the number of RIS elements increases. This suggests that a larger
number of reflecting elements provides more flexibility for the joint phase
shift and UAV placement optimization, leading to better communication quality
and higher gains.

Figure 4: Energy efficiency for different number of RIS elements when $K=4$. TABLE I: Simulation parameters parameter | value | parameter | value
---|---|---|---
$Z_{U}$ | 70 m | $Z_{R}$ | 40 m
$W_{U}^{(0)}$ | (200 m, 50 m) | $W_{R}$ | (200 m, 0 m)
$\lambda$ | 10 cm | $P_{max}$ | 1 W
$d_{r}$ | 5 cm | $d_{c}$ | 5 cm
$B$ | 20 MHz | $\beta_{0}$ | $10^{-2}$
$\alpha_{k}^{UG}$ | 3 | $\alpha_{k}^{RG}$ | 2.4
$\kappa_{k}^{UG}$ | 2 | $\kappa_{k}^{RG}$ | 2
$p_{k}^{c}$ | 1 mW | $p^{r}$ | 1 mW
$m$ | 2 kg | $g$ | 9.8 $\frac{m}{s^{2}}$
$r_{p}$ | 0.2 | $n_{p}$ | 4
$\rho$ | 1.225 $\frac{kg}{m^{3}}$ | $R_{k,min}$ | 100 $\frac{bits}{s}$

##  V Conclusion

In this paper, we proposed an energy-efficient communication scheme for RIS-
assisted UAV networks. By optimizing the On-Off strategy, phase shift of RIS
elements, and the UAV placement and beamforming, our proposed scheme maximized
the energy efficiency of the network. We formulated the optimization problem
as a mixed integer nonlinear program and utilized the BCD approach to tackle
it. We employed a combination of continuous and binary genetic algorithms, as
well as the Adam optimizer, to optimize the various variables of the problem.
Our simulation results showed that the proposed scheme outperforms existing
literature in terms of energy efficiency. Specifically, our scheme achieves
significant improvements in the energy efficiency of downlink communications
between the aerial base station and the GUs.

## References

  * [1] Pogaku, Arjun Chakravarthi, et al. ”UAV-assisted RIS for future wireless communications: A Survey on optimization and performance analysis.” _IEEE Access_ 10 (2022): 16320-16336.
  * [2] M. Javad-Kalbasi and S. Valaee, ”Re-configuration of UAV Relays in 6G Networks,” _2021 IEEE International Conference on Communications Workshops (ICC Workshops),_ Montreal, QC, Canada, 2021, pp. 1-6.
  * [3] X. Cao, B. Yang, C. Huang, C. Yuen, M. Di Renzo, D. Niyato, and Z. Han, ”Reconfigurable intelligent surface-assisted aerial-terrestrial communications via multi-task learning,” _IEEE Journal on Selected Areas in Communications_ , vol. 39, no. 10, pp. 3035–3050, Jun. 2021.
  * [4] S. Li, B. Duo, M. Di Renzo, M. Tao, and X. Yuan, ”Robust secure UAV communications with the aid of reconfigurable intelligent surfaces,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 10, pp. 6402–6417, Apr. 2021.
  * [5] L. Ge, P. Dong, H. Zhang, J.-B. Wang, and X. You, ”Joint beamforming and trajectory optimization for intelligent reflecting surfaces-assisted UAV communications,” _IEEE Access_ , vol. 8, pp. 78 702–78 712, Apr. 2020.
  * [6] X. Liu, Y. Liu, and Y. Chen, ”Machine learning empowered trajectory and passive beamforming design in UAV-RIS wireless networks,” _IEEE Journal on Selected Areas in Communications_ , vol. 39, no. 7, pp. 2042–2055, Jul. 2020.
  * [7] B. Shang, R. Shafin, and L. Liu, ”UAV swarm-enabled aerial reconfigurable intelligent surface (SARIS),” _IEEE Wireless Communications_ , vol. 28, no. 5, pp. 156–163, Oct. 2021.
  * [8] B. Shang, H. V. Poor, and L. Liu, ”Aerial reconfigurable intelligent surfaces meet mobile edge computing,” _IEEE Wireless Communications (Early Access)_ , pp. 1–9, May. 2022.
  * [9] Y. Li, C. Yin, T. Do-Duy, A. Masaracchia, and T. Q. Duong, ”Aerial reconfigurable intelligent surface-enabled URLLC UAV systems,” _IEEE Access_ , vol. 9, pp. 140 248–140 257, Oct. 2021.
  * [10] A. Khalili, E. M. Monfared, S. Zargari, M. R. Javan, N. M. Yamchi, and E. A. Jorswieck, ”Resource management for transmit power minimization in UAV-assisted RIS HetNets supported by dual connectivity,” _IEEE Transactions on Wireless Communications_ , vol. 21, no. 3, pp. 1806–1822, Aug. 2021.
  * [11] M. Javad-Kalbasi, Z. Naghsh, M. Mehrjoo and S. Valaee, ”A new heuristic algorithm for energy and spectrum efficient user association in 5g heterogeneous networks,” _2020 IEEE 31st Annual International Symposium on Personal, Indoor and Mobile Radio Communications,_ London, UK, 2020, pp. 1-7.
  * [12] M. Javad-Kalbasi and S. Valaee, ”Energy and spectrum efficient user association for backhaul load balancing in small cell networks,” _GLOBECOM 2020 - 2020 IEEE Global Communications Conference,_ Taipei, Taiwan, 2020, pp. 1-6,
  * [13] M. Javad-Kalbasi and S. Valaee, ”Centralized and distributed algorithms for energy and spectrum efficient user association in small cell networks,” in _IEEE Transactions on Green Communications and Networking,_ vol. 5, no. 4, pp. 1781-1790, Dec. 2021
  * [14] J. J. Quispe, T. F. Maciel, Y. C. Silva, and A. Klein, ”Joint beamforming and BS selection for energy-efficient communications via Aerial-RIS,” in _Proc. IEEE Global Communications Conference Workshops (GC Wkshps)_ , Madrid, Spain, Dec. 2021.
  * [15] Aung, Pyae Sone, et al. ”Energy-efficient communication networks via multiple aerial reconfigurable intelligent surfaces: DRL and optimization approach.” arXiv preprint arXiv:2207.03149 (2022).
  * [16] Kingma, Diederik P. and Jimmy Ba, ”Adam: A method for stochastic optimization.” arXiv preprint arXiv:1412.6980 (2014).
  * [17] Wei, Zhiqiang, et al. ”Sum-rate maximization for IRS-assisted UAV OFDMA communication systems.” _IEEE Transactions on Wireless Communications_ 20.4 (2020): 2530-2550.
  * [18] Li, Jingyi and Jiajia Liu, ”Sum rate maximization via reconfigurable intelligent surface in UAV communication: Phase shift and trajectory optimization.” _2020 IEEE/CIC International Conference on Communications in China (ICCC). IEEE, 2020._
  * [19] Al-Abiad, Mohammed S. and Md Jahangir Hossain, ”Coordinated scheduling and decentralized federated learning using conflict clustering graphs in fog-assisted IoD Networks.” _IEEE Transactions on Vehicular Technology_ (2022).
  * [20] Gen, Mitsuo and Runwei Cheng, ”Genetic algorithms and engineering optimization,” Vol. 7. John Wiley and Sons, 1999.


# The velocity dependence of dry sliding friction at the nano-scale

Rasoul Kheiri1 Alexey A Tsukanov1,2 and Nikolai V Brilliantov1 1 Skolkovo
Institute of Science and Technology, Moscow, 121205, Russia  2 Research and
Development Centre, TerraVox Global Ltd., Paphos, 8310, Cyprus
<EMAIL_ADDRESS>
###### Abstract

We performed molecular dynamics (MD) experiments to explore dry sliding
friction at the nanoscale. We used the setup comprised of a spherical particle
built up of 32,000 aluminium atoms, resting on a semi-space with a free
surface, modelled by a stack of merged graphene layers. We utilized LAMMPS
with the COMB3 many-body potentials for the inter-atomic interactions and
Langevin thermostat which kept the system at $300K$. We varied the normal load
on the particle and applied different tangential force, which caused the
particle sliding. Based on the simulation data, we demonstrate that the
friction force $F_{\rm fr}$ linearly depends on the sliding velocity $v$, that
is, $F_{\rm fr}=-\gamma v$, where $\gamma$ is the friction coefficient. The
observed dependence is in a sharp contrast with the macroscopic Amontons-
Coulomb laws, which predict the velocity independence of sliding friction. We
explain such a dependence by surface fluctuations of the thermal origin, which
give rise to surface corrugation hindering sliding motion. This mechanism is
similar to that of the viscous friction force exerted on a body moving in
viscous fluid.

##  1 Introduction

According to Amontons-Coulomb laws of friction, for a macroscopic solid object
moving on a surface, dry sliding friction neither depends on the contact are
of the bodies, nor on the sliding velocity. Instead, it is proportional to the
normal load, e.g. [1, 2]. In contrast, the viscous friction force, which
emerges when a body moves through viscous fluid, do depends on the body’s
velocity. The form of the dependence may vary, but normally, it increases with
the increasing speed. Often a linear dependence of the friction force $F_{\rm
fr}$ on the velocity $v$ is observed, $F_{\rm fr}=-\gamma v$, where $\gamma$
is the friction coefficient. Such an important difference stems from different
mechanisms which give rise to the friction force. In the former case, dry
friction originates due to mesoscopic asperities on the surface of macroscopic
bodies; the asperities engage with each other and hinder the motion. In the
latter case, viscous friction emerges due to thermal fluctuation of the local
molecular stress and hence it is determined by the temperature of the system.
Obviously the mechanism of dry friction for macroscopic bodies is athermal,
that is, it is not related to molecular fluctuations [3].

Molecular fluctuations may become important, however, at the meso and
nanoscale, providing a significant impact on the dry friction. At the
mesoscale, studies of the static and dynamic friction and the transition to
sliding demonstrated the importance of elastic deformations, see e.g. [4, 5,
6, 7, 8, 9]; moreover, it has been shown that the onset of sliding may vary
[8]. In Ref. [7] deviations from the Amontons-Coulomb laws have been also
reported – the dependence of friction force on the contact area, and of
friction coefficient on the normal load.

Similarly, recent studies in nanotribology [10, 11, 12] demonstrated that
solid friction at the nanoscale can noticeably deviate from that predicted by
Amontons-Coulomb macroscopic laws. For instance, in the absence of a net
normal load in thin films, friction could be more like the dissipation force
in liquids [13], that is, it can presumable be velocity dependent.

The experimental investigation of dry solid friction at the nanoscale started,
essentially, with the invention of fine experimental tools like atomic force
microscopy (AFM) [14, 15], surface force apparatus (SFA) [14], and quartz
crystal microbalance (QCM) [16, 17, 18] (see also [19, 20]). Still the most
detailed view can be provided by computer experiments – the numerical
investigation on the atomic scale by means of molecular dynamics (MD)
simulation [21]. Since the number of particles simulated in MD is not
”macroscopically” large, one needs to link the molecular system to a somewhat
artificial system – thermostat, which keeps the desired conditions [21]. The
most popular are Anderson [22], Nose-Hoover [23, 24, 25, 26], and Langevin
[27] thermostats. Each of them have some advantages and disadvantages [21];
sometimes thermostating may result in non-physical effects on the interface,
so that further correction are to be applied in the simulations [12, 28, 29].
Hence reliable data may be obtained only under a careful design of the model
[30].

In the present study we analyse dry sliding friction at the nanoscale by means
of MD. We investigate the friction force acting on a spherical nanoparticle
sliding on a solid surface - the boundary of the semi-space, which is modelled
by a several layers of graphene stacked on each other. Molecular fluctuations
of the thermal origin of the free surface yield surface corrugation. This
hinders the motion of the particle on the surface and thus gives rise to the
friction force. Since such a mechanism, associated with the thermal
fluctuations is similar to the friction mechanism of the viscous friction, one
expects that the friction force would be proportional to the sliding velocity
$v$, that is $F_{\rm fr}=-\gamma v$. Where the friction coefficient $\gamma$
may depend on the normal load $F_{\rm N}$.

The rest of the article is organized as follows in the next Sec. II we present
the detail of the MD simulations and the way of obtaining the friction
coefficient from the rough data. In Sec. III we discuss the MD results, with
the focus on the dependence of the friction force on the sliding velocity.
Finally, in Sec. IV we summarize our findings.

##  2 Molecular dynamics simulations

###  2.1 Simulation setup

MD simulations for a spherical particle moving on a fluctuating surface have
been performed using LAMMPS [31], see Figs. 1-2 which provide the sketch of
the model and simulation images. The spherical particle has been made up of
32,000 aluminium atoms, interacting with a classical many-body potential (see
the discussion below). To prepare the aluminum particle we use a cubic FCC
aluminum lattice and then cut from it a sphere.

Figure 1: Sketch of the simulation setup. The tangential, $F_{x}$, and normal
force, $F_{z}$, are exerted on the spherical aluminum particle built up of
32,000 atoms (blue). The semi-space with a free surface is modelled by a stack
of graphene layers (green and yellow), with the fixed bottom layer (yellow).

(a) At the onset of simulation

(b) After some time steps.

Figure 2: Snapshots of the simulated system – aluminium spherical particle
(blue) is moving on several graphene layers (green) to the right. The exerted
normal load is $F_{z}=33.28\times 10^{-9}N$, and the tangential force is
$F_{x}=4.1\times 10^{-9}N$. The temperature is kept around $T=300K$; the
Langvian thermostat is applied.

The semi-space with the free surface has been constructed using a stack of ten
graphene layers merged together. That is, the carbon atoms of the sheets
experience both the inlayer and interlayer interactions, see Fig. 2. Note,
that since each carbon atom participates in four covalent bonds, the
interaction between the surface atoms and the aluminum atoms is relatively
weak. To prevent the vertical motion of the stack of the layers, subject to
the normal load, we fixed the bottom layer (shown yellow in Fig. 1). Moreover,
to prevent the horizontal motion of the layers, subject to tangential stress
we also fixed the side edges of the graphene layers.

For the potential, we use charge-optimized many-body (COMB3) potential. COMB
refers to the second generation of the potential [32], and COMB3 corresponds
to the third-generation of COMB potential [33]. The four terms of the
potential depend on two variables – the charge $q$ and inter-site distance
$r$:

$U[{q},{r}]=U^{es}[{q},{r}]+U^{short}[{q},{r}]+U^{vdW}[{r}]+U^{corr}[{r}],$
(1)

where $U^{es}$ is electrostatic term, $U^{short}$ is charge-dependent short-
range interaction, $U^{vdW}$ is van der Waals interaction, and $U^{corr}$ is
correction term.111For more information see the comb and comb3 pair_styles in
the LAMMPS documents [34] The equilibrium charge on each atom is calculated by
the electronegativity equalization (QEq) method [35]. This potential is of
high accuracy although computationally time consuming.

The COMB3 potential has been utilized to model the interface of graphene and
metals like C-Cu, and C-Al. In the particular case of Al-C, we can refer to
the reference [36] for the graphene-Al interface. In this potential,
interactions between Al and graphene are mostly weekly bonded. That is, the
interface between Al and graphene remains physisorbing. Moreover, we performed
our simulation for $T=300K$, which implies the lack of aluminum
carbide.222Aluminum carbide formation might happen under vacancy defects in
graphene and high enough temperatures. However for the temperatures range of
$300-900K$, without any defects in graphene, the simulations remain without
the formation of aluminum carbide [36]. For simplicity, we set here zero
charges for Al and C atoms; we expect that this would not have any noticeable
impact on the surface fluctuations, which hinder the tangential sliding.

In our simulations we use the Langevin thermostat along with NVE time
integration scheme (Langevin + NVE). This choice has been motivated by the
fact that utilizing the Langevin thermostat supplemented by the NVE
integration, results in the simulations of a canonical ensemble with nearly
conserved energy [27].333It is important to mention that during thermostating,
even in the Langevin dynamics, the removal of the additional energy from the
system is not quite steady. Therefore, an unphysical dissipation might be
introduced in the system [12]. We believe that this kind of setup would be the
best one for an adequate modelling of the surface thermal fluctuations. It
should be noted however, that the properties of surface fluctuations would
depend on the damping parameter of the Langevin thermostat, which is set with
some uncertainty. Indeed, the damping parameter is related to the viscous
properties of the solid material – its bulk dissipative constants, and hence
may be defined from these quantities. Unfortunately, neither the bulk
dissipative constants, for the material built up of graphene sheets, nor the
microscopic theory for such a relation are currently available. Therefore, we
use here plausible values for the damping parameter of the Langevin
thermostat, which implies that the simulation results provide qualitative
behavior of the system. During the current simulations we used the damping
parameter $damp=0.1$.444For more information see the fix Langevin command in
the LAMMPS documents [37]

The simulations begin with the application of the normal load $F_{\rm
N}=F_{z}$ distributed over all aluminium atoms of the sphere, Fig. 1. In
addition the particle gets an acceleration subject to the applied tangential
force $F_{x}$. The equations of motion for all particles are solved by the
Verlet algorithm [21] for a constant temperature of $T=300K$. After some time,
the acceleration is expected to cease and the particle reaches a constant
velocity regime. Here the friction force is equal to the tangential force,
$F_{\rm fr}=F_{x}$. In practice, the relaxation time towards the steady state
may be rather long, therefore we apply some tricks to estimate the steady-
state sliding velocity using relatively short runs; these are discussed below.

In applying different normal and tangential forces we observe that efficient
and adequate simulations require some limitation for the range of the exerted
forces $F_{x}$ and $F_{z}$. To remain in a realistic simulation time, around a
few nanoseconds, we have to choose $F_{x}$ not too large (to reduce the
relaxation time to the steady state) and $F_{z}$ not too small (to increase
the friction force, which again results in the decreasing relaxation time).
Nevertheless, $F_{z}$ cannot be too large, as a very large force would mash
the particle. All in all, we have chosen the intervals, $F_{x}\in[2,5]nN$ and
$F_{z}\in[15,60]nN$.

We simulate a particle moving in the positive $x$-direction, where the
periodic boundary condition have been applied. This allows long-time
simulations. We use the simulation time step of the order of femtosecond
($dt\sim 0.001\times 10^{-12}s$), and the total simulation time was about a
few nanoseconds ($\sim 10^{-9}s$).

###  2.2 Estimation of the steady state velocity

We conjecture that the friction force is proportional to the sliding velocity,
$F_{\rm fr}=-\gamma v$. Consider the equation of motion for a particle moving
on the surface, subject to the tangential force $F_{x}$, friction force
$-\gamma v$ and a random force $\zeta(t)$ with zero mean,
$\left<\zeta(t)\right>=0$,

$m\dot{v}=-\gamma v+F_{x}+\zeta(t).$

After averaging the above equation we obtain the time evolution of the average
velocity, when the initial average velocity is zero, $\langle v(0)\rangle=0$:

$\langle v(t)\rangle=\frac{F_{x}}{\gamma}\left[1-e^{-(\gamma/m)t}\right],$ (2)

which describes the relaxation of the average velocity to its steady value

$\left.\langle v\rangle\right|_{t\to\infty}=v_{st}=\frac{F_{x}}{\gamma}.$

The last equation shows that the friction coefficient may be found from the
(known) tangential force $F_{x}$ and the steady velocity.

Hence, to measure $\gamma$ in computer simulations, one needs to wait until
the steady state is achieved. In reality, however this is not very practical,
since the relaxation time may be rather long. There is a couple of ways to
shorten the simulation time. Firstly, one can use Eq. (2) and apply a fitting
procedure, to find the best fit for $\gamma$ to satisfy this relation.
Secondly, one can utilize a more simple approach, which we use here. Suppose
that we require the accuracy of $\alpha$ (say $\alpha\approx 0.99$) for the
estimate of $\gamma$. Then at the first step we roughly estimate $\gamma$ and
at the second step we predict the duration of the simulation time $t_{s}$,
such that the average velocity $\langle v(t_{s})\rangle$ at $t=t_{s}$
approximates the steady state velocity $v_{st}$, that is, $\langle
v(t_{s})\rangle\approx v_{st}$, with the accuracy $\alpha$. It immediately
follows from Eq. (2) that

$t_{s}\simeq-\frac{m}{\gamma}\,\ln(1-\alpha).$ (3)

In our MD experiments we use the above estimate for the total simulation time.

##  3 Results and discussion

Fig. 2(a) depicts the system at the onset of the simulation, while Fig. 2(b),
after some simulation time, when the surface thermal fluctuations have
developed. During the simulations, the temperature was kept around $300K$.

Fig. 3 illustrates four realizations of the velocity relaxation to a nearly
steady state regime for a constant normal load of $F_{z}=33.28\times 10^{-9}N$
and different horizontal forces $F_{x}$. As it may be seen from the figure,
the total simulation time was about $2ns$, which suffices to obtain good
estimates for the friction coefficient.

Figure 3: Simulated steady-states for the normal load of $F_{z}=33.28\times
10^{-9}N$, and four different horizontal forces $F_{x}$.

Fig. 4 demonstrates the linear dependence of the friction force on the
velocity of the nanoparticle in the steady state sliding. The best fit yields,
$F_{x}=-\gamma v_{st}$, where $v_{st}$ has been obtained as discussed above.

Now we shortly discuss the possible errors in our simulations. As it has been
stated above, one source of errors is the lack of true steady states, and the
approximation of the steady state velocity by its approximate transient value.
The fixation of the edges of the graphene layers violates the symmetry of
surface fluctuations with respect to $x$ and $y$ directions. This may result
in some unphysical effects with an impact on the friction force. Finally, the
Langevin thermostat, as a stochastic setup also contributes to the errors. In
the latter case, however, the error may be mitigated by multiple runs. To
check the reliability of the measured friction coefficient $\gamma$ we perform
a couple of independent runs with the same system parameters. We observe that
the estimate is reliable for the most of the normal forces $F_{z}$, except for
the smallest loads. Here the linearity of friction force on sliding velocity
could be questionable.

Figure 4: Illustration of the linear dependence of friction force on the
steady state sliding velocity for different normal load.

##  4 Conclusion

We investigate numerically, by the means of MD simulation, dry sliding
friction at the nanoscale. In our study we utilize LAMMPS with the most modern
COMB3 many-body potentials for the inter-atomic interactions. In our
simulation setup we use a spherical particles comprised of 32,000 aluminum
atoms, which rests on a semi-infinite space with a free surface. The latter
part of the system is modelled by a stack of merged graphene layers, with a
fixed bottom layer. We use the Langevin thermostat which keep the temperature
at $300K$. We vary the normal load on the particle $F_{z}$ and apply different
tangential force $F_{x}$. We observe that after a transient time, the system
relaxes to a state with a steady sliding velocity $v_{st}$. The steady
velocity is determined by the applied tangential force or, the other way
around, the friction force is determined by the sliding velocity. Hence, based
on our data of the MD experiments we demonstrate that the friction force
linearly depends on the sliding velocity, that is, $F_{\rm fr}=-\gamma v$.
Such a linear dependence holds true for the studied interval of normal and
tangential forces, $F_{z}\in[15,60]\,\,nN$ and $F_{x}\in[2,5]\,\,nN$. It,
however differs from the non-linear dependence expected from the Prandtl-
Tomlinson (PT)model [38]. Although the agreement of our simulation data with
the predictions of the PT model was not expected, the relatively large sliding
velocities used in our study may be among the reasons of the difference.

The main conclusion of our study is that, in contrast to the macroscopic dry
sliding friction, the friction force at the nanoscale is velocity-dependent.
Moreover for the range of normal loads and sliding velocities explored here,
the friction force linearly increases with the increasing velocity. The
physical mechanism of the observed friction behavior is the corrugation of the
contact plane by surface fluctuations of the thermal origin. Such fluctuations
hinder the relative motion of the surfaces, yielding the friction force,
proportional to the sliding velocity. This is similar to the viscous friction
force that acts on a solid body moving in viscous fluid.

AAT and NVB gratefully acknowledge the Russian Foundation for Basic Research
(RFBR), Grant No. 18-29-19198. The research was carried out using the
equipment of the shared research facilities of HPC computing resources at
Lomonosov Moscow State University [39] and Skoltech supercomputer Zhores [40].

## References

## References

  * [1] Popov V L 2010 Contact mechanics and friction (Springer)
  * [2] Berman A, Drummond C and Israelachvili J 1998 Tribology letters 4 95
  * [3] Brilliantov N V, Budkov Y A and Seidel C 2016 Philosophical Transactions A 374 20160143
  * [4] Rubinstein S M, Cohen G and Fineberg J 2004 Nature 430 1005–1009
  * [5] Rubinstein S M, Cohen G and Fineberg J 2006 Physical review letters 96 256103
  * [6] Rubinstein S, Cohen G and Fineberg J 2007 Physical review letters 98 226103
  * [7] Ben-David O, Cohen G and Fineberg J 2010 Science 330 211–214
  * [8] Ben-David O and Fineberg J 2011 Physical review letters 106 254301
  * [9] Lorenz B and Persson B 2012 Journal of Physics: Condensed Matter 24 225008
  * [10] Urbakh M, Klafter J, Gourdon D and Israelachvili J 2004 Nature 430 525–528
  * [11] Mo Y, Turner K T and Szlufarska I 2009 Nature 457 1116–1119
  * [12] Vanossi A, Manini N, Urbakh M, Zapperi S and Tosatti E 2013 Reviews of Modern Physics 85 529
  * [13] Krim J 2012 Advances in Physics 61 155–323
  * [14] Binnig G, Quate C F and Gerber C 1986 Physical review letters 56 930
  * [15] Carpick R W and Salmeron M 1997 Chemical reviews 97 1163–1194
  * [16] Krim J and Widom A 1988 Physical Review B 38 12184
  * [17] Krim J, Solina D and Chiarello R 1991 Physical Review Letters 66 181
  * [18] Krim J 2007 Nano Today 2 38–43
  * [19] Persson B N 2013 Sliding friction: physical principles and applications (Springer Science & Business Media)
  * [20] Mate C M and Carpick R W 2019 Tribology on the small scale: a modern textbook on friction, lubrication, and wear (Oxford University Press, USA)
  * [21] Frenkel D and Smit B 2001 Understanding molecular simulation: from algorithms to applications vol 1 (Elsevier)
  * [22] Andersen H C 1980 The Journal of chemical physics 72 2384–2393
  * [23] Nosé S 1984 The Journal of chemical physics 81 511–519
  * [24] Nosé S 1984 Molecular physics 52 255–268
  * [25] Hoover W G 1985 Physical review A 31 1695
  * [26] Hoover W G 1986 Physical Review A 34 2499
  * [27] Schneider T and Stoll E 1978 Physical Review B 17 1302
  * [28] Benassi A, Vanossi A, Santoro G E and Tosatti E 2012 Tribology Letters 48 41–49
  * [29] Kantorovich L 2008 Physical Review B 78 094304
  * [30] Dong Y, Li Q and Martini A 2013 Journal of Vacuum Science & Technology A: Vacuum, Surfaces, and Films 31 030801
  * [31] Thompson A P and Plimpton S J 2016 Lammps: A general open-source framework for particle-based simulation of materials on multiple scales. Tech. rep. Sandia National Lab.(SNL-NM), Albuquerque, NM (United States)
  * [32] Shan T R, Devine B D, Kemper T W, Sinnott S B, Phillpot S R et al. 2010 Physical Review B 81 125328
  * [33] Liang T, Shan T R, Cheng Y T, Devine B D, Noordhoek M, Li Y, Lu Z, Phillpot S R and Sinnott S B 2013 Materials Science and Engineering: R: Reports 74 255–279
  * [34] LAMMPS 2003-2022 pair_style comb3 command, & pair_style comb3 command URL https://docs.lammps.org/pair_comb.html
  * [35] Rick S W, Stuart S J and Berne B J 1994 The Journal of chemical physics 101 6141–6156
  * [36] Zhang D, Fonseca A F, Liang T, Phillpot S R and Sinnott S B 2019 Physical Review Materials 3 114002
  * [37] LAMMPS 2003-2022 fix langevin command URL https://docs.lammps.org/fix_langevin.html
  * [38] Vanossi A, Manini N, Urbakh M, Zapperi S and Tosatti E 2013 Rev. Mod. Phys. 85(2) 529–552 URL https://link.aps.org/doi/10.1103/RevModPhys.85.529
  * [39] Adinets A, Bryzgalov P A, Voevodin V V, Zhumatii S A, Nikitenko D A and Stefanov K S 2012 Numerical Methods and Programming (Vychislitel’nye Metody i Programmirovanie) 13 160–166
  * [40] Zacharov I, Arslanov R, Gunin M, Stefonishin D, Bykov A, Pavlov S, Panarin O, Maliutin A, Rykovanov S and Fedorov M 2019 Open Engineering 9 512–520