\section{Introduction}
\subsection{Background}
The recent futuristic notion of Metaverse is an extension, a complete simulation, and a mirror of the real world, which is empowered by the maturing high-performance Extended Reality (XR) technologies. Among them, Virtual Reality (VR) provides users with a fully immersive and digital world, which has been used in many fields such as entertainment, socialization, and industry~\cite{metaverse, mobileVR}. To alleviate the obtrusive sense of restriction for users when moving around wearing VR devices, the mobility of these devices is essential. Feasible solutions are applied by using wireless connections and decreasing the weights at the sacrifice of local computation abilities. Thus, even state-of-art VR devices (e.g., HTC Vive~\cite{VRsurvey}) do not have sufficient local computing power to support high-resolution and high-frame-rate applications. In this event, transferring some VR frame generations to a remote server is necessary. Furthermore, under the context of the Metaverse, the boundary between virtual and physical environments becomes more and more blurred, which brings more frequent demands of accessing the virtual world from users. As the Metaverse emphasizes building the bridge between the virtual and real worlds and supporting user inputs anywhere and anytime, it is necessary to consider how to allocate the limited resources to a wide variety of VR users (VUs) with distinct characteristics and demands.

\subsection{Challenges and motivations}
We will first explain the challenges with motivations for this work from the following aspects.

\textbf{User-centric features and user-diverse problems.}
The main reason for its popularity can be contributed to the user-centric services, which serve as the fundamental basis of the next-generation Internet~\cite{metaverse}. The Metaverse is a perceived virtual universe and ideally, allows people to access it anywhere and anytime, with any purpose. The user-centric features also pose a user-diverse problem. The wide variety of use cases (applications) and user-inherent characteristics (e.g., device types, battery power) induce highly different demands for frames per second (FPS), resolutions, etc. Therefore, the first and foremost challenge is how to achieve efficient utilization of those network resources in the multi-user scenario where each user has a different purpose of use and requirements. This propels us to seek a more user-centric and oriented solution to handle widely different users.

\textbf{Wireless communication for VR.} VR is a key feature of an immersive Metaverse socialization experience. Compared to traditional two-dimensional images, generating $360^{\circ}$ panoramic images for the VR experience is computationally intensive. However, as the mobility of VR devices is attached with great importance under the Metaverse context, manufacturers have to lessen the weight at the cost of local computation capability. As a consequence, the existing VR devices lack the local computation ability of high-resolution and frame-rate applications. A feasible solution to powering an immersive socialization experience on VR devices is to make the Metaverse Server help the application frame generation, and send generated frames  to VUs. However, with the frequent demand and high congestion degrees for network resources, it is necessary to lighten the network burden by allocating some devices with higher local computing power to do local generation sometimes. Therefore, this paper considers two cases: (i) server generation, and (ii) local generation.

\textbf{Joint optimization in wireless communication.} In many scenarios, there are more than one important objectives that need to be optimized. In the wireless scenario where two cases of server and local generation are taken into account, two parts are detrimental to data transfer efficiency: (i) The channel access arrangement for VUs (including being assigned with no channel and doing local generation), and (ii) transmission power allocation for VUs. Based on this, this paper focuses on jointly optimizing these variables in the wireless downlink transmission scenario, and considering the VUs' diverse characteristics.

\subsection{Related work and our novelty}
The related work is separated into multiple aspects according to the challenges and motivations since we design our work and make contributions from these multiple domains. The references and our novelties compared to them are expounded in the following.

\textbf{User-centric Metaverse.}
For the time being, people are already very aware of the user-centric particularity in the Metaverse. Lee \emph{et al.}~\cite{metaverse} claim that the Metaverse is user-centric by design, and it will rely on pervasive network access. Consequently, users with diverse purposes, devices, and demands can access the universe anytime and anywhere. Du \emph{et al.}~\cite{usercentricdu} emphasized the user-centric demands in the Metaverse and proposed an attention-aware network resource allocation considering the diversities among users' interests and applications. Both papers above came up with attractive and fresh concepts and discussed the potential challenges and future directions. However, none of them has ever studied a specific scenario with problems and designed novel algorithms to tackle them. In this paper, we design a concrete user-diverse problem scenario and invent a user-centric DRL structure correspondingly.

\textbf{VR over wireless communication.} In recent years, VR services over wireless communication have been thoroughly studied in many previous works. Yang \emph{et al.}~\cite{pengyang} investigated the problem of providing ultra-reliable and power-efficient VR strategies for wireless mobile VUs by Deep Reinforcement Learning (DRL) algorithms. In order to fit the discrete action space required in their DRL algorithm, they quantize the continuous actions into discrete actions. Xiao \emph{et al.}~\cite{chenxuVR} studied the predictive VR video delivering by optimizing the bitrate with DRL methods. Some other works has also demonstrated the excellent performance of DRL methods over wireless communications as its ability to explore and exploit in self-defined environments~\cite{chenxuwireless, chenxuwireless2, FLRL}. However, none of the previous works considered the varying purpose of use and requirements between different VUs, and there are no existing works that have designed a \textit{user-centric} and \textit{user-oriented} DRL method as our proposed solutions.

\textbf{Joint optimization in wireless communications with DRL methods.}
Some remarkable works have investigated using DRL methods to solve joint optimization problems~\cite{MALS}. For instance, Guo \emph{et al.}~\cite{JO1} solved the handover control and power allocation joint problem using Multi-agent Proximal Policy Optimization (MAPPO) and obtained satisfactory results. As mixed continuous-discrete actions can lead to problems of directly embedding the MAPPO algorithm, they use discrete power requirement factors instead of a continuous power allocation to simplify the problem. Thus, the problem they aimed to address does not have heterogeneous actions (both discrete and continuous actions), hence they can directly use MAPPO. He \emph{et al.}~\cite{JO2} researched the joint optimization problems on channel access and power resource arrangements by using DRL methods to determine the channel assignment and conducting traditional optimization methods for power allocation under the channel information. However, none of these works considers a user-diverse scenario or problems with heterogeneous actions. Compared to them, we propose a novel Multi-Agent Deep Reinforcement Learning (MADRL) structure, which is equipped with a user-centric view and able to handle interactive and heterogeneous actions.

\subsection{Methodology and Contributions}
This paper proposes a novel multi-user VR model in a downlink Non-Orthogonal Multiple Access (NOMA) system. Specifically, we optimize the channel access arrangement and downlink power allocation jointly, taking the diversities between VUs into consideration. We designed a novel MADRL algorithm, User-centric Critic with Heterogenous Actors (UCHA), which considers the varying purpose of use and requirements of the users and handles the heterogeneous action spaces. In terms of the backbone of our algorithm, we re-design the widely-used Proximal Policy Optimization (PPO) algorithm~\cite{PPO} with a reward decomposition structure in Critic, and two assymetric Actors.

Our contributions are as follows:

\begin{itemize}
\item \emph{Formulating user-centric VR in the Metaverse:} We study the user-centric Metaverse over the wireless network, designing a multi-user VR scenario where a Metaverse Server assists users in generating reality-assisted virtual environments.
\item \emph{Heterogeneous Actors for inseparable optimization variables:} We created two asymmetric Actors interacting with each other to handle the inseparable and discrete-continuous mixed optimization objectives. Specifically, Actor one is for the channel access arrangement, and Actor two is for the power allocation based on the solution by Actor one.
\item \emph{User-centric Critic for user-diverse scenario:} We crafted a novel user-centric Critic equipped with a more user-specific architecture, in which we decompose the reward into multiple VUs and evaluate the value for each VU. To the best of our knowledge, we are the first to embed the hybrid reward architecture in multi-agent reinforcement learning, and use it to solve communication problems.
\item \emph{Novel and comprehensive simulation:} We conduct extensive experiments and design novel metrics to evaluate our proposed solution. And the experiment results demonstrate the superior performances of UCHA.
\end{itemize}

\subsection{Organization}

The rest of the paper is organized as follows. Section~\ref{sec:models} introduces our system model. Then Sections~\ref{sec:RLenv} and \ref{sec:algorithms} propose our deep reinforcement learning setting and algorithm. In Section~\ref{sec:experiment}, extensive experiments are performed, and various methods are compared to show the prowess of our strategy. Section~\ref{conclude} concludes the paper.

A 6-page short version is under review for the 2023 IEEE International Conference on Communications (ICC)~\cite{ICC}. In that conference version, we demonstrate that the user-centric structure in Critic remarkably accelerates the convergence speed, and using this structure in PPO is much superior to Hybrid Reward DQN~\cite{HRA}. But in that version, we do not consider anything about the resolution, also the optimization of the downlink transmission power, so it is nothing to do with the heterogeneous actions or the MADRL. However, these are all the highlights of this journal version. Besides, this paper also designs new algorithms, metrics, and evaluation methods, compared to the conference version.

\section{System Model}
\vspace{-0.1cm}
\label{sec:models}

Consider a multi-user wireless downlink transmission in an indoor environment, in which $T$ frames are generated by the Metaverse server and sent to $N$ VR Users (VUs) in one second. To ensure a smooth experience, we apply the clock signal from the server for synchronization, and use a slotted time structure, where one second is divided into $T$ time slots (steps), and the duration of each slot is $\iota=\frac{1}{T}$. In each slot, a high-resolution 3D scene frame for each user is generated and sent to $N$ VUs $\mathcal{N}=\{1,2,...,N\}$ via a set of channels $\mathcal{M}=\{1,2,...,M\}$. These VUs have distinct characteristics (e.g., local computation capability) and different requirements (e.g., Frames Per Second (FPS)). This FPS requirement is decided by the specific applications. Assume that one VU is playing VR video games, while another VU is having a virtual meeting. The FPS requirement for the former should be evidently much higher than the latter. Each user can accept VR frame rates as low as a minimum tolerable frame per second (FPS) $\tau_{n,F}$, which is the number of successfully received frames in one second (i.e., $T$ time slots). Our objective is to obtain channel access and downlink power arrangements for each VU.

\subsection{Channel allocation and frame generation}

In terms of channel allocation, we define an $N \times T$ matrix $\boldsymbol{Z}$ such that the element in its $n$th row and $t$th column is $z_n^t$, for $n \in \{1,2,\ldots,N\}$ and $t \in \{1,2,\ldots,T\}$. In other words, $\boldsymbol{Z}$ denotes the selection of downlink channel arrangement. Specifically, in the studied system of one Metaverse server and $N$ VUs, our channel allocation for each VU $n \in \{1,2\ldots,N\}$ includes two cases:
\begin{itemize}
\item \textbf{Case 1: Server-generated frame.} If VU $n$ is assigned for a channel $m$ at time step $t$ (i.e., $z_n^t=m, m\in\mathcal{M}$), the Metaverse server selects the $m$th channel for downlink communication with VU $n$ to deliver the frame that the server generates for VU $n$. In this case, if the sum delay (for generation and transmission) of each frame exceeds slot duration $\iota$, this frame is deemed as a failure.
\item \textbf{Case 2: VU-generated frame.} If VU $n$ is assigned with no channel (i.e.,$z_n^t=0$), it generates the frame locally with a lower computation capability, and at the expense of energy consumption without communicating with the server. In this case, each VU assigned to do local generation will generate the frame of the highest in-time-processable resolution (i.e., can be generated in $\iota$). If the generated resolution is below the minimum acceptable resolution, this frame is deemed as a failure.
\end{itemize}

Therefore, the channel allocation in this paper also includes the decisions on whether the frames for the VUs are generated by the server or the VUs (for simplicity, we sometimes just say channel allocation without explicitly mentioning decisions of \mbox{frame-generation} locations since the former includes the latter). In other words, $z_n^t=m,m\in\mathcal{M}$ indicates VU $n$ is arranged to channel $m$ at time step $t$, and $z_n^t=0$ means user $n$ needs to generate locally. Next, we explain the two cases in detail.

\begin{figure*}[t]
\centering
\setlength{\abovecaptionskip}{-0.4cm}
\includegraphics[width=1\linewidth]{systemmodel.pdf}
\caption{System model. This figure illustrates a single time slot execution, where the frame generation contains two cases, server-generation and local-generation. Then, different metrics are evaluated and rewards are given accordingly. }\vspace{-0.8cm}
\label{fig:model}
\end{figure*}

\subsection{The case of the frame being generated and sent by the server}\label{model:off-Model}
In each time step, the server will manage the downlink channels $\mathcal{M}$ of all VUs $\mathcal{N}$, and it will subsequently allocate the downlink transmission powers $\boldsymbol{p}^t$ with this Channel State Information (CSI). Here we define $\boldsymbol{P}:=\{\boldsymbol{p}^1, \boldsymbol{p}^2, \ldots, \boldsymbol{p}^T\}$ as the downlink power, where $\boldsymbol{p}^t:=[p_1^t, p_2^t, \ldots, p_N^t]$, and enforce $\sum_{n\in\mathcal{N}}p_n^t\leq p_{max}$. As the total delay in this case and the resolutions are important metrics, the achievable rate and resolutions are discussed in the following:

\subsubsection{Achievable rate from the server to each VU}

We adopt the Non-Orthogonal Multiple Access (NOMA) system as this work's propagation model. In our NOMA system, several VUs can be multiplexed on one channel by superposition coding, and each VU exploits successive interference cancellation (SIC) at its receiver. The decoding process follows the approach described by Vaezi~\textit{et~al.}~\cite{vaezi2019multiple}. Specifically, with $\boldsymbol{z}^t$ denoting $[z_1^t,z_2^t,\ldots,z_N^t]$, then among the $N$ VUs, we let $\mathcal{N}_m^t(\boldsymbol{z}^t)$ be the set of $N_m^t(\boldsymbol{z}^t)$ users which receive VR frames from the channel $m$ at time step $t$. Formally, we have
\begin{align}
 \mathcal{N}_m^t(\boldsymbol{z}^t):= \{{n}\in\mathcal{N}~|~{\iffalse \Gamma \fi}z_n^t=m\}, \text{ and } N_m^t(\boldsymbol{z}^t) = |\mathcal{N}_m^t(\boldsymbol{z}^t)|.  \label{eqNm}
\end{align}
The part ``$(\boldsymbol{z}^t)$'' means that $\mathcal{N}_m^t(\boldsymbol{z}^t)$ and $N_m^t(\boldsymbol{z}^t)$ are both functions $\boldsymbol{z}^t$. For simplicity, we sometimes write $\mathcal{N}_m^t(\boldsymbol{z}^t)$ and $N_m^t(\boldsymbol{z}^t)$ as $\mathcal{N}_m^t$ and $N_m^t$. The above practice also applies to other notations throughout this paper.

Based on Page 197 of~\cite{vaezi2019multiple}, we order $N_m^t(\boldsymbol{z}^t)$ VUs of $\mathcal{N}_m^t(\boldsymbol{z}^t)$ such that channel-to-noise ratios between them and the server on channel $m$ are in a decreasing order; formally, recalling that  all $N$ VUs are indexed by $1,2,\ldots,N$, suppose that the above ordering for $N_m^t(\boldsymbol{z}^t)$ VUs of $\mathcal{N}_m^t(\boldsymbol{z}^t)$  produces VU indices $u_1, u_2, \ldots, u_{\mathcal{N}_m^t(\boldsymbol{z}^t)}$; i.e., with $h_{i,{\iffalse \Gamma \fi}m}^t$ denoting the channel attenuation between the server and VU $i$ over channel ${\iffalse \Gamma \fi}m$ at time step $t$, and $(\sigma_{i,{\iffalse \Gamma \fi}m}^t)^2$ denoting the power spectral density of additive white Gaussian noise at VU $i$ over channel ${\iffalse \Gamma \fi}m$ at time step $t$, we have the following\footnote{We write $(\sigma_{i,{\iffalse \Gamma \fi}m}^t)^2$ here for generality. In our experiments of Section~\ref{sec-Experiments}, $(\sigma_{i,{\iffalse \Gamma \fi}m}^t)^2$ is the same $\sigma^2$ for any $i,m,t$.}:
\begin{align}
\begin{array}{l}
\text{define $u_1, u_2, \ldots, u_{\mathcal{N}_m^t(\boldsymbol{z}^t)}$ such that they together form the set $ \mathcal{N}_m^t(\boldsymbol{z}^t)$ of Eq.~(\ref{eqNm})}  \\ \text{\textit{and} }  \frac{|h^t_{u_1,m}|^2}{(\sigma^t_{u_1,m})^2 }\geq \frac{|h^t_{u_2,m}|^2}{(\sigma^t_{u_2,m})^2 }\geq \ldots \geq \frac{|h^t_{u_{\mathcal{N}_m^t(\boldsymbol{z}^t)},m}|^2}{(\sigma^t_{u_{\mathcal{N}_m^t(\boldsymbol{z}^t)},m})^2 } \begin{array}{l} \text{(For the reasoning of such ordering, interested}  \\ \text{\,\,readers can refer to Page 197 of~\cite{vaezi2019multiple}).}   \end{array}
\end{array} \nonumber \\[-20pt]
\label{defineuseq}
\end{align}
Note that if two VUs with indices $u$ and $u'$ result in the same channel-to-noise ratios, their ordering can be arbitrary and this paper places the larger index ahead of the lower index; i.e., $u$ is before (resp., after) $u'$ if $u$ is greater (resp., smaller) than $u'$.  Following the same rationality as\footnote{Page 197 of~\cite{vaezi2019multiple} presents the example of 2 VUs. Our (\ref{defineuseq}) and~(\ref{eq:dlrate}) use the same rationality as~\cite{vaezi2019multiple}. Interested readers can refer to~\cite{vaezi2019multiple} and the references therein for the details.  }~\cite{vaezi2019multiple}, after decoding via SIC, the interferences to VU $n$ satisfying ${\iffalse \Gamma \fi}z_n^t=m$ come from signals sent from server that are intended for users whose positions are before VU $n$ in the sequence $u_1, u_2, \ldots, u_{\mathcal{N}_m^t(\boldsymbol{z}^t)}$; formally, supposing $u_{\nu}$ is $n$ (such $\nu$ exists since $n\in\mathcal{N}_m^t(\boldsymbol{z}^t)$ follows from ${\iffalse \Gamma \fi}z_n^t=m$), those users are $u_1, u_2, \ldots, u_{\nu -1}$ and the sum of the interferences to VU $n$ (i.e., $u_{\nu}$) are given by $\sum_{j=1}^{\nu -1} p_{u_j}^t |h_{n,m}^t|^2$, after we define $p_i^t$ as the transmit power used by server for transmitting the signal intended for VU $i$ at time step $t$.

Then the achievable rate of VU $n$ over its assigned channel ${\iffalse \Gamma \fi}z_n^t=m$  is
\begin{align}
   \hspace{-2pt}{r}_{n}^{t}(\boldsymbol{z}^t,\hspace{-1pt}\boldsymbol{p}^t)\hspace{-1pt}=\hspace{-1pt}W_m\log\hspace{-2pt}\left(\hspace{-2pt}1\hspace{-1pt}+\hspace{-1pt}\frac{{p}_{n}^{t}|h_{n,m}^{t}|^2}{\sum_{j=1}^{\nu -1} p_{u_j}^t |h_{n,m}^t|^2\hspace{-1pt}+\hspace{-1pt}W_m(\sigma_{n,{\iffalse \Gamma \fi}m}^t)^2}\hspace{-2pt}\right)\hspace{-3pt}, \begin{array}{l} \text{for $\nu$ satisfying $u_{\nu}=n$ after defining}\\ \text{related notations via~(\ref{eqNm}) and~(\ref{defineuseq}),}\end{array}\label{eq:dlrate}
\end{align}
where $W_m$ denotes the bandwidth of channel $m$,

\subsubsection{Resolution of server-generated frame}

We denote $D_{n,o}^{t}$ ($n \in \mathcal{N}$) as the transmission data size of the VR frame at time step $t$ that needs to be executed and transmitted by server to user $n$, $\mathcal{G}=\{G_{1}, ..., G_{J}\}$ as the frame sizes with different graphic resolutions in descending order (e.g., $G_1$ is 1440p, and $G_2$ is 1080p), where $G_{1}$ and $G_{J}$ are the lowest and highest acceptable resolutions, respectively. For convenience, we also define $G_{0} = +\infty$ and $G_{J+1}=0$. ${Res}_n^t$ is the received frame resolution of VU $n$ at $t$ (${Res}_n^t \in \mathcal{G}$). Assume that the server-generated frames are all in the highest resolution (which is $G_{1}$), and the transmission data size from server is
\begin{align}
    D_{n}^t = \frac{{Res}_n^t}{{Com}_{n}^{t}}.
\end{align}
In parallel to the fast-developing VR devices, video compression technologies are being developed as well. VR compression leverages the likeness of different images from different cameras and uses advanced slicing and tiling techniques~\cite{compression}. The compression ratio is not always a constant but variable to different VR image qualities and data sizes. To better fit the real situation, we use the $Com_n^t$ as the varying compression ratios of this frame of VU $n$ at time step $t$.

Accordingly, the delay $d_{n}^{t}(\boldsymbol{z}^t,\boldsymbol{p}^t)$ of each frame in time step $t$ is divided into (1) execution time and (2) downlink transmission time:
\begin{align}
    d_{n}^{t}(\boldsymbol{z}^t,\boldsymbol{p}^t) = \frac{D_n^t\times c_n^t}{f_v} + \frac{D_n^t}{{r}_{n}^{t}(\boldsymbol{z}^t,\boldsymbol{p}^t)}, \label{eq:offdelay}
\end{align}
where $f_v$ is the computation capability of the server (i.e., cycles per second), and $c_n^t$ is the required number of cycles per bit of this frame~\cite{COsettings}.

\subsection{The case of the frame being generated by VUs locally}\label{model:local-Model}
When VU is not allocated a channel, it needs to generate the VR frames locally at the expense of \textbf{energy consumption} and \textbf{resolution degeneration} according to its local computing capability (CPU frequency). Let $f_n$ be the computation capability of VU $n$ , and it varies across VUs. Adopting the model from~\cite{energy}, the energy per cycle can be expressed as $e_{n,cyc} = \eta f_n^2$. Therefore, the energy consumption overhead of local computing can be derived as:
\begin{align}
    e_{n,l}^t =
    \begin{cases}
        \mu_n \times D_n^t \times c_n^t \times e_{n,cyc}, &{\iffalse \Gamma \fi}z_n^t \neq 0.\\
        0, &{\iffalse \Gamma \fi}z_n^t =0.
    \end{cases}
     \label{eq:energy}
\end{align}
Inherently, $\mu_n$ is the battery weighting parameter of energy for VU $n$. The battery state of each VU can be different, then, we assume that $\mu_n$ is closer to $0$ with a higher battery.

In terms of the resolution, if VU $n$ is doing local generation, the resolution of the current frame will degenerate to the highest resolution that VU $n$ can process with tolerable delay. Thus, the frame resolution of VU $n$ at $t$ is formulated as:
\begin{align}
    Res_{n}^t =
    \begin{cases}
        G_{1}, &{\iffalse \Gamma \fi}z_n^t \neq 0. \\
        G_{J_{n}^t}, &{\iffalse \Gamma \fi}z_n^t = 0.
    \end{cases}
\end{align}
where
\begin{align}
    &G_{J_{n}^t}/{Com}_n^t \leq f_n \times \iota \leq G_{(J_{n}^t-1)}/{Com}_n^t,~\text{and}~~ 1 \leq J_{n}^t \leq J. \nonumber
\end{align}
Here, $J_{n}^t$ is the rank of the highest available resolution of VU $n$ among all resolution $\mathcal{D}$, and $J$ is the number of resolution types. $\iota$ is the duration of each time step ($\iota = \frac{1}{T}$), and $f_n\times \iota$ is the maximum datasize can be locally generated by VU $n$ in one step.

The overall system model is shown in Fig.~\ref{fig:model}.

\subsection{Problem formulation}
Different users have different purposes for use (video games, group chat, etc.). Therefore, they also have varying expectations of satisfactory FPS $\tau_{n,F}$. For each successive frame, the total delay   the tolerable threshold (occurs in ``Case 1'') and the insufficient resolution ($Res_{n}^t < G_{J}$, occurs in ``Case 2'') lead to a frame failure. We set a frame success flag $I_n^t$ for VU $n$ at step $t$ as:
\begin{align}
    I_n^t =
    \begin{cases}
        0, &\text{if}~z_n^t \in \{1,2,\ldots,M\} \text{ and }d_{n}^{t}(\boldsymbol{z}^t,\boldsymbol{p}^t) > \iota. \\
        0, &\text{if}~z_n^t =0 \text{ and }Res_{n}^{t} < G_{J}. \\
        1, &\text{else}.
    \end{cases}
\end{align}
Our goal is to decide channel allocation and downlink power allocation for the transmission of $T$ frames. The objectives can be divided into 1) fulfilling the FPS requirements of different users as possible and minimizing the total frame transmission failure, 2) optimizing VUs' device energy consumption regarding their battery state, and 3) increasing VUs' received frames resolutions.

We define an indicator function $\chi[A]$ which takes $1$ if event $A$ occurs and takes $0$ otherwise. From the above discussion, our optimization problem is
\begin{align}
\max\limits_{\boldsymbol{Z},\boldsymbol{P}}
&\left \{ \omega_1 \left ( \min\limits_{n\in\mathcal{N}} \left[\left(\sum_{t=1}^T I_n^t\right)-\tau_{n,F}\right] \right) + \frac{1}{N} \sum_{n=1}^N\sum_{t=1}^T[-\omega_2 e_{n,l}^t + \omega_3 {Res}_n^t] \right \} \label{eq:objectivefunction}\\
s.t.~C1:&{\iffalse \Gamma \fi}z_n^t\in\{0,1,...,M\},~\forall{n}\in\mathcal{N},\forall{t}\in\{1,2,\ldots,T\},\\
     C2:&\sum_{n\in \mathcal{N}} \chi[{\iffalse \Gamma \fi}z_n^t \neq 0] p_n^t \leq p_{max}, ~\forall{t}\in\{1,2,\ldots,T\}.
\end{align}

The $\omega_1,\omega_2, \omega_3$ are the weighting parameters of frame failures, local device energy consumption, and VU devices received frames resolution. In practice, these parameters will be reflected in the reward setting~\ref{sec:RLenv}. The first part of the objective function (i.e., $\min\limits_{n\in\mathcal{N}} [(\sum_{t=1}^T I_n^t)-\tau_{n,F}]$) is to make every VU fulfill their FPS requirements as possible, and the second part is to minimize the local energy usage and maximize the frame resolutions. Constraint $C1$ is our integer optimization variable which denotes the computing method and channel assignment for each user at every time step. Constraint $C2$ is the limit of the downlink transmission power for all VUs.

\subsection{The execution flow of the system}

Here we describe the execution flow of the system given a solution to the optimization problem. The key idea of the two cases design (server or local generation) is to utilize the limited network resources and let some VUs spare the resources by doing local generation at each time. Thus, this paper uses a slotted time structure by applying the clock signal from the server for synchronization. Specifically, in each time slot $t$, one frame of each VU will be executed. Considering the very short duration $\iota$ of each time slot, we assume that the channel attenuation $h_{n,m}^t$ varies across different slots but remains the same in one slot. At the beginning of each slot $t$ (frame), the server collects the time-varying information of each VU (e.g., channel attenuation, compression rate), and sends the decisions on channel access (local or server generation, if server generation, which channel) to all VUs through a dedicated channel. Then, VUs assigned to no channel will locally generate the frame immediately, and the others will wait for the server to transmit the frames to them.

\subsection{Motivation of using MADRL to solve the optimization problem}

The formulated problem is non-convex and NP-hard, and the optimization variables: computation method choice, channel arrangement, and power allocation make the formulated problem highly coupled with inseparable mixed-integer non-linear programming (MINLP) optimization problems, where the discrete variable $\boldsymbol{\Gamma}$ and continuous variable $\boldsymbol{P}$ are inseparable. Moreover, this formulated problem is \textbf{time-sequential}, where the number of variables increases with $T$. Thus, the traditional optimization methods are unsuitable for our proposed problem due to the daunting computational complexity. Also, as the problem contains too many random variables, \textit{model-based reinforcement learning (RL)} approaches that require transition probabilities are infeasible techniques to tackle our proposed problem. \textit{Heuristic search} methods can possibly be a solution to sequential problems. However, it doesn't change the approximation of the policy, while naively making improved action selections given the current value function. Besides, a huge tree of possible continuations is usually considered when using this method~\cite {RLintro}, which further makes it impractical to apply heuristic search in such a complicated scenario with a huge dimension of decisions. Therefore, it is highly necessary to design a comprehensive \mbox{model-free} Multi-Agent Deep Reinforcement Learning (MADRL) method to tackle the problem with heterogeneous optimization variables and distinct requirements from all VUs.

\section{Environment in Our Multi-Agent Deep Reinforcement Learning (MADRL)}
\label{sec:RLenv}
To tackle a problem with DRL method, designing a comprehensive reinforcement learning environment based on the formulated problem is the first and foremost step. For a reinforcement learning environment, the most important components are (1) State: the key factors for an agent to make a decision. (2) Action: the operation decided by an agent to interact with the environment. (3) Reward: the feedback for Agent to evaluate the action under this state. Thus, we expound on these three components next.

\subsection{State}
In the DRL environment, weeding out less relevant and less time-varying variables is essential. Therefore, we set two agents: $Agent_1$ for optimizing the channel arrangement, and $Agent_2$ for allocating the downlink transmission power.

\subsubsection{$Agent_1$ State $s^t_1$} We included the following attributes into the state: (1) Each VU's frame size: $D_n^t={Res}_n^t/{Com}_n^t$. (2) Each VU's remaining tolerable frame transmission failure count. (3) The channel attenuation of each VU: $h_{n,m}^t$. (4) The remaining number of time slots: $(T-t)$. (5) Gap from requirement to each VU: $\tau_{n,F}-\sum_{i=1}^tI_n^i$.

\subsubsection{$Agent_2$ State $s^t_2$} Only after obtaining the CSI, the power allocation can be finished based on it. As a result, the \textbf{action} of $Agent_1$ is significant to the decision of $Agent_2$. Besides the action, the other important attributes in $s^t_1$ should be considered by $Agent_2$ as well. Therefore, we use the concatenation as the state: $s_2^t=\{a_1^t;s_1^t\}$.

\subsection{Action}
The appropriate action settings that are directly related to the optimization variables are critical for finding a good solution. In this environment, the action $a_1^t$ and action $a_2^t$ are respectively the $Agent_1$ and $Agent_2$ actions, explained below.  \\
\begin{figure}[t]
    \centering
    \setlength{\abovecaptionskip}{-1cm}
    \includegraphics[width=0.8\linewidth]{actionencode.pdf}
    \caption{UL action encoding method.}
    \label{fig:actionencode}
    \vspace{-0.5cm}
\end{figure}
\subsubsection{$Agent_1$ Action $a_1^t$}
The discrete action of $Agent_1$ is the channel allocation: $\boldsymbol{z}^t = \{z_1^t, z_2^t,..., z_N^t\}$. In a DRL environment, the discrete actions are actually discrete numbers (indicators), then, we need to give consecutive discrete indicators to each action. As $\boldsymbol{z}^t$ contains $N$ elements ($N$ VUs), and each element has $M+1$ possible values ($M$ channels and plus 1 for decision on the cases), we use the $N$-bit \mbox{$(M+1)$}-number code and decode it into the decimal indicator, as shown in Fig.~\ref{fig:actionencode}. The action is given by:
\begin{align}
    a_1^t = \sum_{n=1}^{N}z_n^t(M+1)^{n-1},~\forall t\in \{1,\cdots,T\}.
\end{align}
\subsubsection{$Agent_2$ Action $a_2^t$}
The continuous action of $Agent_2$ should be the downlink transmission power $\{p_1^t,p_2^t,...,p_N^t\}$. However, it is impractical for the RL agent to allocate power for all VU within the sum power constraint. Therefore, we add a softmax layer to the $Agent_2$. Accordingly, the $a_2^t$ turns into the portions of the $p_{max}$. We use $P_n^t$ denotes the portion allocated to $n$ at $t$:
\begin{align}
    a_2^t = \{P_1^t, P_2^t,...,P_N^t\},~\forall t\in \{1,\cdots,T\}.
\end{align}

\subsection{Reward}
In the traditional CTDE framework~\cite{CTDE}, the rewards are shared by different agents. However, in our scenario, the objectives and optimization variables are highly different. e.g., If VU $n$ is allocated to do the local computation, the rewards for energy consumption and resolution degeneration are irrelevant to $Agent_2$. Therefore, we design the rewards for $Agent_1$ and $Agent_2$ separately. Moreover, as we consider a user-centric scenario, we decompose the rewards for $Agent_1$ and $Agent_2$ among different VUs, and construct a novel algorithm in Section~\ref{sec:algorithms}.

\subsubsection{$Agent_1$ Reward for each VU $R_1^t[n]$}
The $Agent_1$ is responsible for channel access arrangement, which takes an important position in the whole system. The rewards for any VU $n$ in $Agent_1$ contain: (1) a descending reward $R_r^t[n]$ for different received frame resolutions from high to low. (2) a penalty for every transmission failure $R_f^t[n]$. (3) a weighted penalty $R_e^t[n]$ for energy consumption corresponding to VU's battery state: $\omega_e\times e_{n,l}^t(\mu_n)$. To fulfill our objective in Eq.~(\ref{eq:objectivefunction}), we give (4) a huge reward/penalty $R_{w}^t[n]$ according to the \textbf{Worst VU} at the final step: $\omega_{end} \times \left ( \min\limits_{n\in\mathcal{N}} \left[\left(\sum_{t=1}^T I_n^t\right)-\tau_{n,F}\right] \right)$. However, the sparse reward is devastating to a goal-conditioned DRL environment without taking any actions~\cite{sparsereward}, and only giving a huge penalty at the final step can make it very hard for training. To avoid the sparse reward and accelerate training, we set an \textbf{early termination flag}. For any VU, if the frame failure times exceed the tolerant failure times, (failure times of $\text{VU}_n$ $> T-\tau_{n,F}$), we assign an (5) additional early termination penalty $R_{term}^t[n]$ according to the number of left frames: $\omega_f \times (T-t)$, and end this episode immediately. $\omega_e, \omega_{end}, \omega_f$ above are both hyperparameters to be set in the experiments.

\subsubsection{$Agent_2$ Reward for each VU $R_2^t[n]$}
The $Agent_2$ is responsible for the downlink power allocation. We weed out some rewards from $R_1^t[n]$ that are not related to $Agent_2$. Therefore, the rewards for every VU $n$ in $Agent_2$ contain: (1) different resolution reward: $R_{r}^t[n]$, (2) transmission failure penalty $R_f^t[n]$ (3) the huge penalty for Worst VU: $R_{w}^t[n]$, and (4) the early termination (task fail) penalty $R_{term}^t[n]$.

For every reward, we narrow the range by dividing the number of VUs $N$ to ease the training.

\section{Our User-centric MADRL Approach} \label{sec:algorithms}
Our proposed User-centric Critic with Heterogeneous Actors (UCHA) structure uses the state-of-the-art Proximal Policy Optimization (PPO) algorithm as the backbone. Inspired by the effective Hybrid Reward Architecture (HRA)~\cite{HRA}, we design a user-centric Critic, which evaluates the current state-value in a more user-specific view. And considering the heterogeneous action space that incorporates both discrete and continuous actions, we create a heterogeneous Actor structure. Thus, the preliminaries, PPO and HRA will first be introduced. We will then explain UCHA.

\subsection{\textbf{Preliminary}}
\subsubsection{Backbone: Proximal Policy Optimization (PPO)}
\textbf{Why PPO?} As we emphasize developing a \textit{user-centric} model which considers VUs' varying purpose of use and requirements and uses multiple agents with heterogeneous actions, the algorithm's policy stability and the ability for dealing with both discrete and continuous actions are essential. Therefore, the Advantage Actor-Critic structure based algorithms~\cite{A2C} that directly evaluate the V values for states instead of Q values for actions is an advisable choice. Among them, the Proximal Policy Optimization (PPO) by openAI~\cite{PPO} is an enhancement of the traditional Advantage Actor-Critic which fulfills the two requirements with better sample efficiency by using two separate policies for sampling and training, and it is more stable by applying the policy constraints.

PPO has been actively used in solving wireless communication problems and Metaverse~\cite{wfiot, MECRLsurvey}, and its prowess has been demonstrated in many scenarios, like the recent widely discussed chatbot, \mbox{ChatGPT}~\cite{ChatGPT}. Next, we will introduce the pipeline of PPO and expound on its two pivotal features: (i) Importance sampling, and (ii) Policy constraint.

\textbf{Use of importance sampling.} Importance sampling (IS) refers to using another distribution to approximate the original distribution~\cite{owen2000safe}. In order to increase the sample efficiency, PPO uses two separate policies (distributions) for training and sampling to better utilize the collected trajectories, which uses the theory of importance sampling.~\cite{PPO}. To distinguish between the two policies, we use $\pi_\theta$, $\pi_{\Bar{\theta}}$ to denote the policies for training and sampling, where $\pi$ is the policy network and $\theta, \Bar{\theta}$ are the parameters. In practice, the Temporal-Difference-1 (TD1) state-action pairs are used, then, the objective function is reformulated as:
\begin{align}
    J(\theta) =\mathbb{E}_{(s^t,a^t)\sim\pi_{\theta}} \left [\pi_\theta(s^t,a^t) A^t \right ]
    =\mathbb{E}_{(s^t,a^t)\sim\pi_{{\Bar{\theta}}}} \left [\frac{\pi_\theta(s^t,a^t)}{\pi_{{\Bar{\theta}}}(s^t,a^t)} A^t \right ]
    \approx \mathbb{E}_{(s^t,a^t)\sim\pi_{{\Bar{\theta}}}} \left [\frac{\pi_\theta(a^t|s^t)}{\pi_{{\Bar{\theta}}}(a^t|s^t)} A^t \right ], \label{eq:IS}
\end{align}
$A^t$ is short for the advantage function $A(s^t,a^t)$ (a standard term in reinforcement learning to measure how much is a certain action a good or bad decision given a certain state). How to sensibly estimate the $A$ has been widely discussed. In this paper, the advantage function is estimated by the truncated version of generalized advantage estimation (GAE)~\cite{GAE}, which will be explained in ~\ref{sec:UCHA}. We assume $\pi_{\theta}(s^t) = \pi_{\Bar{\theta}}(s^t)$ here as calculating the probabilities of states occurrence is impractical~\cite{PPO}, and we use $\approx$ instead of $=$ to avoid misunderstandings.

\textbf{Add KL-divergence penalty.}
To increase the stability, we need to reduce the distance between the two distributions $\theta$ and $\Bar{\theta}$. Therefore, Trust Region Policy Optimization (TRPO)~\cite{TRPO}, which is the predecessor of PPO, uses a Kullback-Leibler (KL) divergence constraint to the objective function to limit the distance between the two distributions by directly setting $D_{KL}(\pi_{\theta}||\pi_{{\Bar{\theta}}}) \leq \varrho$, where $\varrho$ is a hyper-parameter which will be set in their experiments. Nonetheless, this constraint is imposed on every observation, and it is very hard to use in practice. Thus, PPO re-formulates its objective function into the following~\cite{PPO}:
\begin{align}
    \Delta\theta = \mathbb{E}_{(s^t,a^t)\sim\pi_{{\Bar{\theta}}}}[\triangledown f^t(\theta,A^t)],~\text{and}~
    f^t(\theta, A^t)=\min\{\frac{\pi_{\theta}(a^t|s^t)}{\pi_{\Bar{\theta}}(a^t|s^t)}A^t,
    clip(\frac{\pi_{\theta}(a^t|s^t)}{\pi_{\Bar{\theta}}(a^t|s^t)}, 1-\epsilon, 1+\epsilon)A^t\}.\label{eq:actorobj}
\end{align}

\textbf{Critic loss.} Considering the particularity of our Critic structure, we separate the Critic and Actor networks, instead of using the shared layers like the implementation in the famous Stable Baselines3 (SB3)~\cite{SB3}. The loss function is as follows:
\begin{align}
    L(\phi) = (V_\phi(s^t)-V^t_{target})^2,~\text{and}~V^t_{target} = A^{GAE}+\gamma V_{\phi'}(s^{t}). \label{eq:criticloss}
\end{align}
The $V$ is the Critic (value) network, and $\phi$ is the parameters. Then, $V_{\phi}(s)$ means the state-value~\cite{RLintro} of state $s$ generated by Critic network. Note that $\phi'$ is the parameter of the target Critic network, and it will be replaced periodically by $\phi$. This prevailing trick is to increse the stability of the target~\cite{RLintro}. Besides, the advantage function is estimated by the GAE algorithm, and we simply use $A^{GAE}$ here.

\subsubsection{Hybrid Reward Architecture (HRA)}

Since there are always mixed objectives in the communication scenarios that lead to hybrid rewards reflected in reinforcement learning, we usually need to find a solution to tackle the hybrid reward problem. In this scenario, the hybrid rewards refer to the different rewards for different VUs, as we consider their distinct inherent characteristics and requirements. This issue of using RL to solve a high dimensional objective function was first studied in~\cite{HRA}. In their work, they proposed the HRA structure for Deep Q-learning (DQN) which aims to decompose rewards into different reward functions according to each objective. e.g., Decomposing rewards into objectives for energy consumption and delay minimization, and calculating two different losses accordingly. The illustration of HRA is shown in Fig.~\ref{fig:hra}. HRA can exploit domain knowledge to a much greater extent and has remarkable effects on environments with different roles and objectives, which has been demonstrated to speed up learning with better performances in various domains such as video game playing~\cite{HRAgame}, mimicry learning, complementary task~\cite{HRAmimic}, etc. However, to our best known, no one has explored this structure in depth and used it in solving communication problems. This paper re-design the normal Critic into a user-centric Critic by embedding HRA structure. Specifically, we decompose the rewards into different VUs instead of different objectives and calculate the losses of each VU accordingly. Next, we will expound on our novel and effective algorithm, User-centric Critic with Heterogeneous Actors (UCHA).

\begin{figure}[t]
\centering
\includegraphics[width=0.85\linewidth]{HRA.pdf}
\vspace{-1cm}
\caption{The Hybrid Reward Architecture from~\cite{HRA}. The reward is decomposed into multiple domains, and different from the normal value network which only outputs one Q value for the action, HRA outputs multiple Q values for every domain.}\vspace{-1cm}
\label{fig:hra}
\end{figure}

\subsection{User-centric Critic with Heterogeneous Actors (UCHA)}
\label{sec:UCHA}
In contrast to decomposing the overall reward into separate sub-goal rewards as done in HRA, we built a user-centric reward decomposition Critic, which takes in the rewards of different users and calculates the actions-values separately. In other words, we give the network a view of the value for each user, instead of merely evaluating the overall value of an action based on an overall state. Simultaneously, considering the highly different roles, action spaces, and objectives of the two agents, we design the Heterogeneous Actors structure. The two agents are updated by different rewards and advantages, which will be expounded in the following.

\textbf{Function process:} In each episode, when the current transmission is accomplished with the channel access arrangement $a_1^t$ from $Agent_1$, and downlink power allocation $a_2^t$ from $Agent_2$, the environment will issue every VU's rewards $\boldsymbol{R}_1^t=\{R_1^t[1], R_1^t[2],..., R_1^t[N]\}$ for $Agent_1$, and $\boldsymbol{R}_2^t=\{R_2^t[1],R_2^t[2],...,R_2^t[N]\}$ for $Agent_2$ as feedback to different VUs. The global state $s_1^t$ and the next global state $s_1^{t+1}$ will be sent to the user-centric Critic to generate the state values for each VU. Then, the rewards $\boldsymbol{R}_1^t$ and $\boldsymbol{R}_2^t$ with state-values will be used to calculate the advantages $\boldsymbol{A}_1^t=\{A_1^t[1], A_1^t[2],..., A_1^t[N]\}$ and $\boldsymbol{A}_2^t=\{A_2^t[1], A_2^t[2],..., A_2^t[N]\}$ for $Agent_1$ and $Agent_2$. The user-centric Critic takes in the global state and advantages for calculating the Critic losses $\{L_1,L_2,...,L_N\}$for different VU, and updates with the sum losses similar to HRA. Note that the \textbf{global} state ($s_1^t$), reward ($\boldsymbol{R}_1^t$), and advantage are all from $Agent_1$, as those of $Agent_2$ are mainly a portion extracted from $Agent_1$. In terms of the Actors, we use $\boldsymbol{A}_1^t$ and $\boldsymbol{A}_2^t$ to update Actors of $Agent_1$ and $Agent_2$, separately. The above-mentioned process is illustrated in Fig. \ref{fig:alg}.

\begin{figure}[t]
\centering
\includegraphics[width=1\linewidth]{algorithm.pdf}
\vspace{-1cm}
\caption{The structure of User-centric Critic with Heterogeneous Actors (UCHA). The top of the figure is the function process and overall map of UCHA, while the bottoms are the network structures in practice (coding). The two distinct Actors take in different states and select different actions, and the Critic evaluates the global state, outputting values for each VU, separately.}
\label{fig:alg}
\vspace{-0.8cm}
\end{figure}

\textbf{Heterogeneous actors update function:} In equation~(\ref{eq:actorobj}) (i.e., $\Delta\theta = \mathbb{E}_{(s^t,a^t)\sim\pi_{{\Bar{\theta}}}}[\triangledown f^t(\theta,A^t)]$), we established the policy gradient for PPO Actor, and our UCHA  uses the sum-advantages of every VU for the update. Then, we have the gradient $\Delta\theta_1$, $\Delta\theta_2$ for $Agent_1$ and $Agent_2$ as:
\begin{align}
    \Delta\theta_1 = \mathbb{E}_{(s^t,a^t)\sim\pi_{{\theta}'}}[\triangledown f^t(\theta,(\sum_{n=1}^N A_1^t[n])],~ \Delta\theta_2 = \mathbb{E}_{(s^t,a^t)\sim\pi_{{\theta}'}}[\triangledown f^t(\theta,(\sum_{n=1}^N A_2^t[n])].\label{eq:gradient}
\end{align}
where $A_1^t[n], A_2^t[n]$ denote the advantages of different VUs for $Agent_1$ and $Agent_2$. A truncated version of generalized advantage estimation (GAE)~\cite{schulman2015high} is chosen as the advantage function:
\begin{align}
    &A_1^t[n] = \delta_1^t[n] + (\gamma\lambda)\delta_1^{t+1}[n]+...+(\gamma\lambda)^{\bar{T}-1}\delta_1^{t+\bar{T}-1}[n],\\
    &A_2^t[n] = \delta_2^t[n] + (\gamma\lambda)\delta_2^{t+1}[n]+...+
    (\gamma\lambda)^{\bar{T}-1}\delta_2^{t+\bar{T}-1}[n],
\end{align}
where
\begin{align}
    \delta_1^t[n]=R_1^t[n]+\gamma V_{\phi'}(s_1^{t+1})[n]-V_{\phi'}(s_1^t)[n],~\delta_2^t[n]=R_2^t[n]+\gamma V_{\phi'}(s_1^{t+1})[n]-V_{\phi'}(s_1^t)[n].
\end{align}
$\bar{T}$ specifies the length of the given trajectory segment, $\gamma$ specifies the discount factor, and $\lambda$ denotes the GAE parameter. And these parameters will be specified in experiments. $V_{\phi'}$ is the target Critic with parameters $\phi'$, which will be periodically replaced by $\phi$. The user-centric Critic has $N$ outputs, where $V_{\phi'}(\cdot)[n]$ means the Critic-head for VU $n$, which takes in the state, and outputs the state-value for VU $n$. Note that we both use $s_1^t$, but with their separate rewards, because $s_1^t$ is the global state, and we only give the Critic an overall view of the whole task.

\textbf{User-centric Critic update function:} In this user-centric Critic, we compute the value losses for each VU separately to enable  UCHA to have a user-specific view. Similar to the update method in HRA~\cite{HRA} that uses the sum losses of different components, we use the sum losses of different VUs. Therefore, the traditional PPO-Critic loss in Eq.~(\ref{eq:criticloss}) (i.e., $L(\phi) = (V_\phi(s^t)-V^t_{target})^2$) is re-formatted into:
\begin{align}
    L(\phi) =\sum_{n=1}^N L(\phi)[n]= \sum_{n=1}^N \left (V_{\phi}(s_1^t)[n]- V_{target}^t[n]\right )^2, \label{eq:HRloss}
\end{align}
where $V_{target}[n] = A_1^t[n]+\gamma V_{\phi'}(s_1^{t})[n]$, and $L(\phi)[n]$ is the Critic loss for VU $n$.

As explained above, this user-centric Critic uses centralized training with equation~(\ref{eq:HRloss}). The computation complexity with real-time expenditure in experiments will be discussed in Section~\ref{sec:experiment}.

\begin{figure}[!t]
    \vspace{-15pt}
        \renewcommand{\algorithmicrequire}{\textbf{Initiate:}}
        \renewcommand{\algorithmicensure}{\textbf{Output:}}

        \begin{algorithm}[H]
            \caption{\label{alg:PPO} User-centric Critic with Heterogeneous Actors (UCHA)}
            \begin{algorithmic}[1]
                \REQUIRE Actor 1 parameter $\theta_1$, Actor 2 parameter $\theta_2$, Hybrid critic parameter $\phi$ and target network $\phi'$, initial state $s_1^0$.
                \FOR{iteration = $1,2...$}
                    \STATE $Agent_1$ executes action according to $\pi_{\theta_1^{'}}(a_1^t|s_1^t)$
                    \STATE Get reward $R_1^t[1],R_1^t[2],...,R_1^t[N]$ and the current state of $Agent_2$: $s_2^t$
                    \STATE $Agent_2$ executes action according to $\pi_{\theta_2^{'}}(a_2^t|s_2^t)$
                    \STATE Get reward $R_2^t[1],R_2^t[2],...,R_2^t[N]$ and the next state of $Agent_1$: $s_1^{t+1}$, and $s_1^t \leftarrow s_1^{t+1}$
                    \STATE Sample \{$s^t, a^t, (R_1^t[1],R_1^t[2],...,R_1^t[N]), (R_2^t[1],R_2^t[2],...,R_2^t[N]),s^{t+1}$\} till end
                    \STATE Compute advantages \{$A_1^t[1],A_1^t[2],...,A_1^t[N]$\} for $Agent_1$, \{$A_2^t[1],A_2^t[2],...,A_2^t[N]$\} for $Agent_2$, and target values\{$V_{target}^t[1],V_{target}^t[2],...,V_{target}^t[N]$\}
                    \FOR{$k$ = $1,2,...,K$}
                        \STATE Shuffle the data's order, set batch size $bs$
                        \FOR{$j$=$0,1,...,\frac{\textbf{Trajectory length}}{bs}-1$}
                            \STATE Compute gradients for Actor 1 and Actor 2 using Eq.~(\ref{eq:gradient}).
                            \STATE Update Actors separately by gradient ascent
                            \STATE Compute Value losses for each VU
                            \STATE Update Critic with MSE loss using Eq.~(\ref{eq:HRloss})
                        \ENDFOR
                        \STATE Assign target network $\phi' \leftarrow \phi$ every $C$ steps
                    \ENDFOR
                \ENDFOR
            \end{algorithmic}
        \end{algorithm}
        \vspace{-1.5cm}
\end{figure}

\section{Experiments} \label{sec-Experiments}
In this section, we conduct extensive experiments to compare UCHA with the baselines, highlighting the remarkable performance achieved by it. We will first introduce our baseline algorithms and metrics, and then present our results with detailed analysis.
\label{sec:experiment}

\subsection{\textbf{Baselines}}
As we introduce a novel UCHA algorithm with unique actor and critic structures, we weed out the key points of UCHA one by one, and design the following baseline algorithms:
\begin{itemize}
    \item \textbf{Independent PPO (IPPO)}. The most intuitive way of using RL in such a cooperative interactive environment is to implement two independent RL agents interacting with each other. We implemented two independent PPO with user-centric Critics for optimizing the channel access and downlink power allocation. They both use a normal and separate Actor and Critic and update with their own states, actions, and rewards. \textbf{This baseline is to examine the effect of the heterogeneous Actors structure in UCHA}.
    \item \textbf{Heterogeneous-Actor PPO (HAPPO)}. We also implement an HAPPO structure that is similar to UCHA. The HAPPO does not apply the user-centric architecture, but merely uses one normal Critic that outputs a single value for the whole global state. Therefore, in HAPPO, we use the sum rewards of each Agent as the single reward. i.e., $Agent_1$ uses rewards $\boldsymbol{R_1^t}=\sum_{n=1}^N R_1^t[n]$, and $Agent_2$ uses $\boldsymbol{R_2^t}=\sum_{n=1}^N R_2^t[n]$. \textbf{This baseline is to testify the effect of the User-centric structure of the Critic}.
    \item \textbf{Random}. The two random Agents select actions randomly, which represents the system performance if no optimization strategy is performed. \textbf{The random policy serves as a naive baseline to show the results if no optimization has been conducted.}
\end{itemize}
\vspace{-0.2cm}

\subsection{Metrics}
We introduce a set of metrics to evaluate the effectiveness of our proposed methods.

\begin{itemize}
    \item \textbf{Worst VU frames}. We define the Worst VU frames, i.e., $\min\limits_{n\in\mathcal{N}}(\sum_{t=1}^TI_n^t-\tau_{n,F})$ in Eq.(\ref{eq:objectivefunction}), as an important metric to show the distance between the successfully received frames and the required FPS of the worst VU.
    \item \textbf{Successful FPS}. The number of successful frames among total $T$ frames determines the FPS of the VR frames and hence fluidity of the Metaverse VR experience.
    \item \textbf{The sum local energy consumption}. We first illustrate the sum energy consumption during training, and then investigate it for different VUs in detail.
    \item \textbf{The received frame resolutions}. The received frame resolutions of different VUs will be shown to testify the comprehensive view of our proposed methods.

\end{itemize}
\vspace{-0.2cm}

\subsection{Channel Attenuation}

When $z_n^t=m$, the channel attenuation becomes a very important variable for the downlink rate. In this paper, the channel attenuation is simulated as $h_{n,m}^t = \sqrt{\beta_n^t}g_{n,m}^t$. The Rician fading is used as the small-scale fading, $g_{n,m}^t = \sqrt{\frac{K}{K+1}}\bar{g}_{n,m}^t + \sqrt{\frac{1}{K+1}}\tilde{g}_{n,m}^t$, where $\bar{g}_{n,m}^t$ stands for the Line-Of-Sight (LOS) component, and $\tilde{g}_{n,m}^t$ means the Non-LOS (NLOS) component that follows the standard complex normal distribution $\mathcal{C}\mathcal{N}(0,1)$. The large-scale fading $\beta_n^t = \beta_0 (L_n)^{-\alpha}$, and $L_n$ represents the distance between the $n$th VU and server. $\beta_0$ denotes the channel attenuation at the reference distance $L_0 = 1$ m. The path-loss exponent $\alpha$ is simulated as $2$, and the Rician factor $K$ is simulated as $3$. Note that as the duration of each time slot is too short, we don't consider the geographical movements of VUs in each time step.

\subsection{Numerical Setting}
Consider a $30\times30$ $m^2$ indoor space where multiple VUs are distributed uniformly across the space. We set the number of channels to be $3$ in each experiment configuration, and the number of VUs to be from $5$ to $8$ across the different experiment configurations. To simplify the simulation, we use the most common standard of the resolutions but ignore the
various horizontal FoVs and aspect ratios for different VR devices. We set the resolutions of one frame can be 1440p ($2560\times1440$, which is known as QHD or Quad HD), 1080p ($1920\times1080$, which is known as Full HD), 720p ($1280\times720$, which is known as HD), and below (considered as insufficient resolution). Each pixel is stored in 16 bits~\cite{oculus} and the factor of compression is selected in the uniformly random distribution $[300,600]$~\cite{compression}. And each frame consists of two images for two eyes. Therefore, the data size of each frame can be $D_n^t\in\{\frac{2560\times1440\times16\times2}{compression},\frac{1920\times1080\times16\times2}{compression}, \frac{1280\times720\times16\times2}{compression}\}$. The maximum flashed rate, $T$ frames in one second is taken to be 90, which is considered a comfortable rate for VR~\cite{VRsurvey} applications. The bandwidth of each channel is set to $10\times180$ kHZ. The required successful frame transmission count $\tau_{n,F}$ is uniformly selected from $[60,90]$, which is higher than the acceptable of $60$~\cite{VRsurvey}. The maximum powers of the server are $100$ Watt. The server computation frequency is $10$ GHz and the computation capabilities of each VU are selected uniformly from $[0.3, 0.9]$ GHz. For all experiments, we use a single NVIDIA GTX 2080 Ti and $2\times 10^5$ training steps, and the evaluation interval is set to be $500$ training steps. As there are several random variables in our environment (e.g., channel attenuation, compression rates, user location), all experiments are conducted under \textbf{ten different global random seeds}, and the error bands are drawn to better illustrate the model performances.

\subsection{Train-time result analysis}
Reasonable and comprehensive experiments are of vital importance in such a complicated scenario. In this section, we will illustrate and analyze some evaluated metrics completely during training. For brevity, we show the performances of different models in two experimental configurations: the less complicated scenario with 6 VUs and the more complicated one with 8 VUs. And the overall results of every scenario are shown in Table~\ref{table:results}. We also analyze the computational complexity, and then the value losses of each VU are illustrated to give a more ``user-specific'' analysis.
 Note that to better evaluate the different metrics, we do not apply early termination if the task fails in the evaluation stage.

\begin{figure*}[t]
\centering
\subfigtopskip=1pt
\subfigbottomskip=1pt

\subfigure[Training reward with 6 VUs.]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=1\linewidth]{experiment/36reward.pdf}
\label{fig:36reward}
\vspace{-10mm}
\end{minipage}
}%
\subfigure[Energy consumption with 6 VUs.]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=1\linewidth]{experiment/36energy.pdf}
\label{fig:36energy}
\vspace{-10mm}
\end{minipage}%
}%
\subfigure[Worst VU frames with 6 VUs.]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=1\linewidth]{experiment/36worstvu.pdf}
\label{fig:36worstvu}
\vspace{-10mm}
\end{minipage}
}%

\subfigure[Training reward with 8 VUs.]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=1\linewidth]{experiment/38reward.pdf}
\label{fig:38reward}
\vspace{-10mm}
\end{minipage}
}%
\subfigure[Energy consumption with 8 VUs.]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=1\linewidth]{experiment/38energy.pdf}
\label{fig:38energy}
\vspace{-10mm}
\end{minipage}%
}%
\subfigure[Worst VU frames with 8 VUs.]{
\begin{minipage}[t]{0.33\linewidth}
\centering
\includegraphics[width=1\linewidth]{experiment/38worstvu.pdf}
\label{fig:38worstvu}
\vspace{-10mm}
\end{minipage}
}%

\caption{Train-time model performances with 6 and 8 VUs (3 channels). In the simpler configuration (6 VUs), UCHA has a much quicker convergence rate than other baselines, and in the more complicated configuration (8 VUs), UCHA obtains much more remarkable performance in all aspects. }\vspace{-1cm}
\label{fig:train}
\end{figure*}

\subsubsection{Train-time model performance in ``$3\text{ Channels, }6\text{ VUs}$'', ``$3\text{ Channels, }8\text{ VUs}$'' configurations} The training reward, sum local energy consumption, and worst VU frames show an overall upward trend as training progresses. When pitted against these metrics, UCHA performed the best out of the tested baseline algorithms in both configurations.

In the experimental setting with $6$ VUs, although HAPPO (without the user-centric structure in Critic) is able to attain similar peak rewards in the end training stages, UCHA converges in about one-fifth of training steps taken for HAPPO to achieve convergence, as shown in Fig.~\ref{fig:36worstvu}. As for the energy usage shown in Fig.~\ref{fig:36energy}, both UCHA and HAPPO increase the energy consumption (i.e., do more local generation) first, and decrease it in the following steps, while the training rewards both increase steadily. Apparently, both UCHA and HAPPO learned to decrease the congestion degree in each channel by allocating some VUs for local generation successfully. They initially tried to decrease congestion degrees to fulfill the requirements of the FPS by allocating more VUs for local generation, and then saving the local battery energy by decreasing the local generation times gradually. There is a sudden drop in the training reward of UCHA at about $100$k steps (Fig.~\ref{fig:36reward}), which is assumed to be the sharp decrease in energy usage (local computation times). This means that UCHA tried to further save energy usage by a fast decline in the computation times, but too many VUs sharing the limited channel and downlink power resources can be devastating to the overall transmission process. On balance, UCHA and HAPPO both succeed in lowering local energy usage to a similar level as if no optimizations are conducted (random policy). As a matter of fact, lowering energy usage to this random policy level is impressive. Because intuitively, we need to sacrifice more local energy to fulfill the task, and the reward and worst VU frames they obtained are extremely higher than the random policy.

Different from UCHA and HAPPO, IPPO fails to find an acceptable solution even in this simpler scenario, but we can still observe from Fig.~\ref{fig:36reward} that the reward is increasing during training. IPPO has the worst stability across different global random seeds, which is reflected in the huge error bands and the sudden drops even in the late training stages. The above observations signify that the Agents in IPPO fail to work cooperatively and they don't achieve good arrangements on channel access and downlink power compared to UCHA and HAPPO. However, IPPO still has a good performance in the simplest $3\text{ Channels, }5\text{ VUs}$ scenario. (Table~\ref{table:results})

The results of the experiment with $8$ VUs show that UCHA is extremely superior to HAPPO and IPPO in almost all aspects. Compared to the simpler configuration with $6$ VUs, UCHA used many more steps to reach the peak reward, as shown in Fig.~\ref{fig:38reward}. Moreover, in this scenario, the local energy consumption (local generation times) is increasing all the time for UCHA, which is different from the observations in the $3\text{ Channels, }6\text{ VUs}$ configuration. It is reasonable because the scenario with more VUs needs more local computation arrangements to avoid unacceptable interference. And the rationality of doing so is reflected in the steadily increasing reward (Fig.~\ref{fig:38reward}) and the slower rise of the energy usage (Fig.~\ref{fig:38energy}). The complete results are shown in Table~\ref{table:results}. We can observe that UCHA obtains the best performance for almost every metric under every scenario. This demonstrates that decomposing the reward and using sum-losses which provides a user-centric view for the RL agent, is a good approach to tackling the formulated problem.
\begin{table}[t]
\vspace{-10pt}
\centering
\caption{Overall results\vspace{-10pt}}
\label{table:parameter}
\vspace{-0.2mm}
\scalebox{0.6}{
    \begin{tabular}{cccccc}
        \hline
        Scenario & \makecell{Energy \\usage (in Joules)} & \makecell{Worst VU:\\$\min(\sum\limits_{t=1}^T I_n^t-\tau_{n,F})$} & Reward & \makecell{Train step\\time (ms)} & \makecell{Execution step\\time (ms)}\\ \hline
        \multicolumn{6}{c}{UCHA} \\ \hline
        $3\text{ Channels, }5\text{ VUs}$ & $11.78$ & $9.97$ & $80.04$ & $19$ & $0.95$\\
        $3\text{ Channels, }6\text{ VUs}$ & $10.84$ & $4.76$ & $69.72$ & $18$ & $0.94$\\
        $3\text{ Channels, }7\text{ VUs}$ & $23.67$ & $5.43$ & $62.78$ & $20$ & $0.94$\\
        $3\text{ Channels, }8\text{ VUs}$ & $37.25$ & $2.56$ & $65.6$ & $22$ & $0.98$\\ \hline
        \multicolumn{6}{c}{HAPPO} \\ \hline
        $3\text{ Channels, }5\text{ VUs}$ & $14.46$ & $4.33$ & $75.46$ & $15$ & $0.92$\\
        $3\text{ Channels, }6\text{ VUs}$ & $13.37$ & $-1.53$ & $64.98$ & $15$ & $0.91$\\
        $3\text{ Channels, }7\text{ VUs}$ & $9.88$  & $-43.76$ & $16.10$ & $16$ & $0.90$\\
        $3\text{ Channels, }8\text{ VUs}$ & $21.96$ & $-60.67$ & $-2.83$ & $18$ & $0.96$\\ \hline
        \multicolumn{6}{c}{IPPO} \\ \hline
        $3\text{ Channels, }5\text{ VUs}$ & $10.76$ & $1.32$ & $72.22$ & $20$ & $0.55$\\
        $3\text{ Channels, }6\text{ VUs}$ & $9.41$ & $-49.67$ & $14.46$ & $27$ & $0.53$\\
        $3\text{ Channels, }7\text{ VUs}$ & $10.96$ & $-66.34$ & $-16.81$ & $32$ & $0.55$\\
        $3\text{ Channels, }8\text{ VUs}$ & $12.23$ & $-73.83$ & $-14.00$ & $56$ & $0.57$\\ \hline
    \end{tabular}
}
\label{table:results}
\vspace{-1cm}
\end{table}

Then we illustrate the value losses for each VU (calculated by the user-centric Critic) to show the convergence of our proposed UCHA in Fig.~\ref{fig:vloss}. Overall, all losses of VUs are decreasing, which shows the convergence of the hybrid Critic even with eight branches (one branch for one VU). However, the range and the decline patterns of the losses are all different. Here, we list four representative VUs for detailed discussion. In this setting, $\text{VU}_1$ is the toughest one for optimizing as it has no sufficient local computing power for a high resolution (i.e., higher than 720p), which is doomed to frame loss, and it has a high FPS requirement (80 FPS). Therefore, $\text{VU}_1$ must be arranged carefully, or it can cause huge losses. $\text{VU}_2$ has a high loss initially, as it has a high requirement of FPS. However, this can be easily tackled by allocating it with more local computation as $\text{VU}_2$ has sufficient local computing power for 1440p resolution. $\text{VU}_6$ is in a more involved situation. It can generate 1080p frames locally and has low battery energy. Thus, if it is arranged to do local computing, although there is no frame loss, it will be assigned energy consumption and resolution degeneration penalties. On the contrary, $\text{VU}_3$ can generate 1440p frames locally with high frame failure tolerance, and its loss is always low.

\begin{figure}[t]
\centering
\vspace{-0.5cm}
\includegraphics[width=0.9\linewidth]{experiment/vloss.pdf}\vspace{-0.4cm}
\caption{The Critic losses (Value losses) for each VU when using UCHA in 3 channels 8 VUs scenario.}\vspace{-0.8cm}
\label{fig:vloss}
\end{figure}

\textbf{Computational complexity.} We  analyze the computational complexity and show the exact time by the figure. We use $K^l_{A}, K^l_{A'}, K^l_{C}$ to denote the number of neurons in layer $l$ of Actor1, Actor2, and the user-centric Critic. And $(\mathcal{A}_0, {\mathcal{A}'_0}, {\mathcal{C}_0})$, $(L_{A}, L_{A'}, L_C)$ be the size of the input layers (proportional to the state dimension shown in Fig.~\ref{fig:alg}), and the number of training layers of the three parts. Each training step contains two Actor training and one Critic training, and considering the mini-batch size $B$ in the training stage, we have the complexity in one training step as{$O(B(\mathcal{A}_0K^1_{A}+\sum_{l=1}^{L_{A}-1}K^l_{A}K^{l+1}_{A} + \mathcal{A}'_0K^1_{A'}+\sum_{l=1}^{L_{A'}-1}K^l_{A'}K^{l+1}_{A'} + \mathcal{C}_0K^1_{C}+\sum_{l=1}^{L_{C}-1}K^l_{C}K^{l+1}_{C}))$. And according to~\cite{complexity}, the computational complexity depends on the total number of convergence steps to the optimal policy. In practice, the network training can be performed offline for a finite number of steps at a centralized-powerful unit (such as the server). We use Table~\ref{table:results} as an intuitive illustration for the time of a single training and execution (in milliseconds) step.

\subsection{Evaluating performance of the user-specific effect}
In this section, we load the trained models of UCHA, HAPPO, and IPPO, evaluate them in the same evaluation environment. Noted that to better compare different algorithms,  we weed out the early termination constraint (i.e., if any VU's frame loss times reach the tolerant limit, this episode ends immediately) in evaluation, which means all evaluations finish the whole $90$ frames transmissions.

\begin{figure*}[t]
\centering
\subfigtopskip=1pt
\subfigbottomskip=1pt
\vspace{-0.5cm}
\subfigure[The received frame resolutions of each VU. The table on the left shows the VU inherent characteristics (i.e., Max resolution when local computing, Battery state, and FPS demand), and the right bar chart compares the received resolutions frames of each VU in one second when using different algorithms.]{
\begin{minipage}[t]{0.9\linewidth}
\centering
\includegraphics[width=1\linewidth]{experiment/resolutions.pdf}
\label{fig:resolution}
\vspace{-5mm}
\end{minipage}
}%

\subfigure[The local computation times and energy usage. The color of the VU index denotes the battery states.]{
\begin{minipage}[t]{0.9\linewidth}
\centering
\includegraphics[width=1\linewidth]{experiment/userenergy.pdf}
\label{fig:userenergy}
\vspace{-10mm}
\end{minipage}%
}%
\vspace{-0.25cm}
\caption{Detailed user-specific evaluation.}\vspace{-0.8cm}
\label{fig:train}
\end{figure*}

Fig.~\ref{fig:resolution} presents the frame resolutions of each VU, in which the left table is the inherent characteristics of each VU (i.e., Max resolution when local computing, battery state, and FPS requirement), and the right bar chart compares the resolutions of each VU with different algorithms. And Fig.~\ref{fig:userenergy} further shows the local computation times and local energy usage of each VU. To simplify the illustration, we only set battery states as Low, Middle, and High, and a fixed local max resolution for each VU in this evaluation. Noted that the ``fail'' in the bar chart means this frame is allocated to be processed on the server, but not transmitted to VU in time, and the resolution less than 720p is also deemed as a failure.

To begin with, only UCHA is able to fulfill the FPS requirements of different VUs, and it does obtain a ``user-specific'' ability. For those with less local computing power (e.g., $\text{VU}_1$, $\text{VU}_2$, $\text{VU}_6$) and with low batter energy (e.g.,$\text{VU}_5$), UCHA avoid allocating them to do many local computations. We can observe from Fig.~\ref{fig:userenergy} that UCHA decides to allocate $\text{VU}_1$ and $\text{VU}_5$ to do remote computing every time, as $\text{VU}_1$ with a high FPS requirement doesn't have the sufficient local computing power for 720p (i.e., the lowest resolution deemed as successful), and $\text{VU}_5$ has a low battery and a high local max resolution. Therefore, $\text{VU}_5$ can consume more energy with the high local resolution if it is allocated to computing locally. On the contrary, $\text{VU}_3$,$\text{VU}_7$,$\text{VU}_8$ are arranged to do more local computation by UCHA, as they have high battery energy and good local computation capabilities, which can spare the channel and downlink power resources to other VUs. Furthermore, we notice that $\text{VU}_4$ receives the least frames successfully, but this is reasonable, as $\text{VU}_4$ is in a less FPS-aware scenario and has the lowest FPS requirement, and at some steps, we need to allocate more resources to other VUs with higher requirements. The above remarkable performances demonstrate that UCHA has a ``user-centric'' ability.

Different from UCHA, HAPPO and IPPO fail to do well in such a user-centric task, and IPPO performs worse than HAPPO. We can observe from the local computation arrangement in Fig.~\ref{fig:userenergy} that HAPPO (without the user-centric Critic) learns a very extreme policy. $\text{VU}_2$, $\text{VU}_3$, and $\text{VU}_7$ are allocated to do local computation almost every time, as $\text{VU}_3$ and $\text{VU}_7$ have high battery energy which leads to relatively low local computing penalty. Although $\text{VU}_2$ has a low battery, it only generates 720p resolution locally, and the energy usage is lower. In terms of IPPO, the disorderly and unreasonable allocations signify that the agents in it are non-cooperative and do not achieve an overall good channel and power selection.\vspace{-5pt}

\section{Conclusion\vspace{-5pt}}
\label{conclude}
In this paper, we study a user-centric multi-user VR for the Metaverse over wireless networks. Users with varying requirements are considered, and a novel user-centric DRL algorithm called UCHA is designed to tackle the studied problem. Extensive experimental results show that our UCHA has the quickest convergence speed and achieves the highest reward than other algorithms, and UCHA has successfully gained the user-specific view. We envision our work to motivate more research on calibrating deep reinforcement learning for research problems of the Metaverse.

{
\begin{spacing}{1.16}
\renewcommand{\refname}{~\\[-25pt]References\vspace{-5pt} }


\section{Introduction}

In the race to reduce power consumption and increase processing capability, the use of spin rather than charge promises a new generation of microelectronics. Spintronics based on metallic multilayers employs these structures for data storage and spin transport for information exchange. Conventional spintronics employs spin currents and spin densities that are due to non-equilibrium electrons that are still close to the Fermi energy.~\cite{Hirohata_2014}

Even though metal-based spintronics is a very well established field, it continues to rapidly evolve. Among the current goals of spintronics research are the manipulation of spin ensembles carried by electrons that are energetically farther away from the Fermi energy.
For non-equilibrium electrons in general, scattering events determine the electronic transport in magnetic materials, regardless of the source that drives the currents, for instance, electrical spin injection from magnetic (metallic or semiconductor) electrodes \cite{Okamoto2016,Tombros2007,johnson2015}, or manipulating spin polarization by alternating-current (AC) magnetic fields via Zeeman interaction \cite{Zhang2018,Busl2010}, to name but a few.
If one considers electron injection at high energies or, in particular, optical excitation of electronics, dynamics with an essential electronic energy dependence and its interplay with spin transport becomes more important. For this hot-electron spin transport an intermediate ``superdiffusive'' regime was identified~\cite{PhysRevLett.105.027203}, which has characteristics somewhat different from ballistic and diffusive transport. Based on experimental evidence~\cite{Rudolf2012, PhysRevLett.117.147203}, it is believed that superdiffusive spin currents can be launched from a ferromagnetic layer into adjacent metallic layers and can contribute to the transfer of spin angular momentum that is needed in the ultrafast demagnetization of ferromagnets. An alternative method to describe superdiffusive spin currents, is provided by a particle in-cell approach~\cite{PhysRevB.98.224416,PC_binder}. This method to solve the spin-dependent Boltzmann equation can be relatively easily adapted to ab-initio input.

In this work, we focus on the dynamics and transport of electrons created far away from the Fermi energy by fs laser irradiation. We study a prototypical Fe/Au bilayer system, where optical pulses, which drive the hot-electron spin currents, are absorbed in both the ferromagnetic and normal metal layers. Our approach is based on the Monte Carlo model for spin-dependent electron dynamics developed in Ref.~\cite{Briones_2022}. In addition to the effects of different collision processes, such as secondary electron generation and elastic scattering, we investigate the influence of different thicknesses of the magnetic layer and the interface transmissions with the intent to understand and optimize spin currents in the structure.

\section{Theoretical Approach}

The aim of our study is to predict how different interactions influence the nonequilibrium dynamics and spin injection in FM/M bilayer systems. In our approach, we consider free electron states above $E_\mathrm{F}$ as essentially free and focus on the influence of high-energy electrons in spin transport. This section presents the algorithms, equations and set of parameters that will be used for the simulations presented in this manuscript.

\begin{figure}[H]

    \includegraphics[scale =0.7, center]{graphs2/FeAu_scheme.eps}%

    \caption{Schematic representation of the laser irradiation and scattering events in the Fe/Au bilayer system. There are two possible interactions to consider: electron–atom interactions where only the direction of flight changes (elastic scattering) and electron-electron interaction generating secondary electrons (impact ionization). When electrons enter in contact with the interface, their spin determines their transmission probability, and they cross accordingly.}
    \label{fig:Scheme}

\end{figure}

\subsection{Monte Carlo method}\label{Sec::MC_method}

The asymptotic Monte Carlo trajectory method \cite{MC_asymp} is a statistical technique that models binary collision interactions by sampling several trajectories until an estimate of the possible outcomes is obtained. Probability theory is used to implement an algorithm for random sampling of variable $x$. In probability theory, all possible events $p(x)$ are integrated into a variable referred to as the cumulative distribution function (CDF) $F(x)$ in order to construct a formula that can be used to obtain a value for the variable $x$ \cite{MCparticle}. A more thorough and accurate description of how the Monte Carlo method is done can be found in literature \cite{alireza, Landau_MC,HUTHMACHER2015242}. Here, we will study the dynamics of excited particles using the same technique explained in our previous work \cite{Briones_2022}. During the simulation, any electron interaction process is treated by random sampling. The time between two successive collisions (time of free flight) $\tau$ can be sampled with the random variable $\mathrm{R} \in [0,1]$ as
\begin{equation}\label{eq:freetime}
 \tau = -\nu^{-1}_0\log (\mathrm{R}),
\end{equation}
when the scattering rate $\nu_0$ is constant. However, assuming a constant total scattering rate, independent of energy, is a statistical overestimation. To compensate, we introduce a further possibility with an energy-dependent probability, which allows the particle to continue its trajectory unperturbed. For the case of several scattering mechanisms as studied in this work, we perform a random sample of collisions using the probability function $p(x)$ which can be replaced with either the differential cross-section $d\sigma/d\Omega$, scattering rates $\nu$ or characteristic times $\tau$.

\subsection{Material parameters}\label{subsection::Laser_parameter}

We consider a Fe/Au bilayer material irradiated by a Gaussian-like laser pulse of $6\;\mathrm{eV}$ photon energy. The laser irradiates the material from the iron side, as shown in Fig. \ref{fig:Scheme} and the maximum laser pulse irradiation is centered at $\SI{0}{\femto\second}$. The iron layer thickness is finite, but the gold layer thickness, is taken to be infinite and the electrons are tracked only in the first $\SI{10}{\nano\meter}$. Excited electrons then interact within the material, either changing only their momentum (elastic scattering) and/or generating secondary electrons (impact ionization). As soon as an electron reaches the interface, its transmission to the other layer is determined by its spin and generally also its energy.

There are two major parameters for this simulation that we tune and study: (i) The thickness of the Fe layer,
and (ii) the transmission probability $T$ at the interface. We focus on the excitation of high-energy electrons where the transmission probability is rather insensitive to the electronic energy
(see figure 4 in \cite{PhysRevLett.119.017202}). Thus, we approximate the transmission probability as a constant.
and also assume it to be independent of propagation direction of the electrons reaching the interface
i.e., it is assumed to be identical for
electrons traveling from iron to gold and the other way around. The constant transmission probability  will be taken taken from Ref. \cite{PhysRevLett.119.017202} as $T^{\uparrow} = 0.95$, and $T^{\downarrow} = 0.25$ for spin-up and spin-down electrons, respectively.

\section{Scattering processes}

We trace the dynamics of optically excited electrons by modeling pure jump processes~\cite{FoundationMC, SimulationMC,HUTHMACHER2015242}. We consider two kinds of interactions in this study, namely the electron-atom interaction, which is treated as an elastic scattering process here, and the electron-electron impact ionization, in which energy is transferred to a secondary electron through an inelastic process.

Elastic scattering processes will  be modeled as a deflection  without transfer of energy. As indicated in subsection \ref{Sec::MC_method}, the probability function $p(x)$ will be replaced by the differential cross-section based on Mott's cross-section in dependence of the scattering angle. We will use $\tau_{el} = 25\;\mathrm{fs}$ as the elastic scattering time~\cite{PhysRevB.71.233104}, which was applied to study spin transport and spin dynamics in Refs.~\cite{PhysRevB.98.224416,PhysRevB.85.235101}.

Electron impact ionization is a process in which incident electrons ionize electrons from occupied states into unoccupied states. The newly generated high-energy electrons are referred to as secondary electrons, whereas the optically excited electrons are referred to as primary electrons. With the energy lost by the primary electron of $\Delta E$  and the binding energy $I$ of the bound electron, the final energy of the newly ionized electron $E_s$ above Fermi level is $E_s = \Delta E - I$. The amount of transferred energy is assumed to be half of the energy of the incident electron as it was done in Ref.~\cite{RITCHIE19651689}. For the inelastic scattering rate ($\tau_{ee}^{-1}$) we will consider the energy-dependent collision rate of an excited electron at temperature $T_{e}$~\cite{PhysRevB.87.035139}. Both collision processes and their inclusion are described in more detail in Ref. \cite{Briones_2022}.

\section{Results}

The behavior of the nonequilibrium spin transport and spin dynamics will be analyzed in a  Fe/Au bilayer system after fs-laser excitation. The laser pulse irradiates the sample from the iron layer (see Fig. \ref{fig:Scheme}), and the distance to the interface is measured from this point inwards. Therefore, the thickness of the iron layer will determine the location of the boundary between these materials. Throughout this work, we will use a variety of iron thicknesses, and they will be specified when required. As for the thickness of the gold layer, it is always infinity and we trace the dynamics within the first
$\SI{10}{\nano\meter}$.
The effects on the spin transport can be addressed by studying the spin currents generated at the interface of the material. The spin current density $j_s$ is defined as
\begin{equation}\label{eq::Spcrr}
    j_s(z,t)\propto q[ \langle \eta^{\uparrow} v_\uparrow \rangle - \langle \eta^{\downarrow} v_\downarrow \rangle ]

\end{equation}
where $q$ is the charge of the electron, $\eta^{\uparrow}$ ($\eta^{\downarrow}$) and $v_\uparrow$ ($v_\downarrow$) are the particle density and the velocity for spin up (spin down), respectively. An average over the number of electrons at a specific volume and time is denoted by the angle brackets. Both expressions in eq. \eqref{eq::Spcrr} can be separated as:

\begin{eqnarray}
    J_{up} &\propto  &\langle \eta^{\uparrow} v_\uparrow \rangle \label{eq:jup}\,,\\
    J_{down} &\propto  &\langle \eta^{\downarrow}  v_\downarrow \rangle\,\label{eq:jdown}.
\end{eqnarray}

\subsection{Spin current density in bulk iron}

We start by examining the behavior of the spin current density in bulk iron. We then extend the analysis to a more complex system, namely the Fe/Au bilayer. As it was analyzed in Ref.~\cite{Briones_2022}, particles reach an average kinetic energy of less than $\SI{2}{\eV}$ very rapidly because secondary electrons are
likely to be generated during the initial time period of the nonequilibrium dynamics.

\begin{figure}[H]
\centering
    \includegraphics[scale =0.68, center]{graphs2/heatmapDepth_Iron31nm.eps}

    \caption{Development of the spin current density $j_s$, defined in eq. \eqref{eq::Spcrr}, in bulk iron over time and in depth, considering secondary electron generation. The dashed lines indicate the position where the spin current density is evaluated further.}
    \label{fig:Spin_current}

\end{figure}

First, let us focus on the irradiation of bulk iron by femtosecond laser pulses as it was described in reference \cite{Briones_2022}. A simulation is performed using bulk-like Fe film which is excited by a Gaussian-like fs laser pulse with a photon energy of $\SI{6}{\eV}$.

As a result of the laser excitation, majority and minority electrons are brought out of equilibrium with different kinetic energy distributions, and also random initial directions. We use open boundary conditions at the end of the calculation region of $\SI{30}{\nano\meter}$ thickness, allowing electrons to escape completely from the simulated system.

Figure \ref{fig:Spin_current} shows the evolution of the spin current density, as defined by eq.~\eqref{eq::Spcrr}, as a function of time and depth. Note that the magnitude of spin current density remains positive over time and throughout the material. This indicates that, on average, majority electrons move towards the depth of the material. In addition, the magnitude of spin current density increases steeply with increasing depth and reaches a maximum value at a around $\SI{10}{\nano\meter}$, which will be analyzed later.

A more detailed observation of the trend in the spin current density is provided in fig. \ref{fig:etavel_J_BL}. Here, we depict the contributions to the spin current density in time at different depths, as indicated in fig. \ref{fig:Spin_current}. The top graph represents the spin current density for both spin up (lhs) and spin down (rhs) electrons, the middle graph depicts the average velocity $\langle v \rangle$, and the bottom graph displays the average particle density $\langle \eta \rangle$, terms used in eqs. \eqref{eq:jup} and \eqref{eq:jdown}. The color gradient indicates at what depth the variables were selected for analysis. The variables shown in fig. \ref{fig:etavel_J_BL} all seem to exhibit the same behavior in time, but a different maximum value for different depths. We observe again how the average particle density $\langle \eta \rangle$ shows a homogeneous density distribution further in time, which was already discussed in ref. \cite{Briones_2022}.

\begin{figure}[H]
\centering
    \includegraphics[scale =0.7, center]{graphs2/etaveljupdown_IronSEAbsBL.eps}

    \caption{An analysis of the evolution of the spin up (lhs) and spin down (rhs) current density $J$ (top graph), average velocity $v$, and particle density $\eta$ in bulk iron, as defined in eq. \eqref{eq::Spcrr}. The physical quantities are depicted for different depths at the location of the dashed lines in fig. \ref{fig:Spin_current}.}
    \label{fig:etavel_J_BL}

\end{figure}

Based on these simulations, we have identified two possible causes for this behavior, either the type of scattering or the absorption profile. In fig. \ref{fig:MaxJ_iron}, it is shown how the maximum values of spin current density $j_s$ are affected by the type of scattering and absorption profile. The red curve has been calculated using Beer-Lambert's absorption probability to excite primary electrons and including secondary electron generation. The black curve was done by keeping the same absorption profile, but considering only optically excited electrons. The magnitude of the maximum spin current density $j_s$ in both graphs increases with depth. This increment then slows down until a certain depth (around $\SI{10}{\nano\meter}$), after which it gradually decreases. These two graphs then differ only qualitatively, which is in agreement with ref. \cite{Briones_2022}, where it is demonstrated that secondary electrons affect spin current density quantitatively rather than qualitatively. As we consider that the absorption profile influences the observed behavior at different maxima of spin current density, we need to examine the effect of a different profile on primary electrons. This particular scenario involves a homogeneous absorption profile, i.e., there will be the same probability of absorption at every depth.

\begin{figure}[H]
\centering
    \includegraphics[scale =0.68, center]{graphs2/Maxj_thickness_Iron.eps}

    \caption{Maximum values of the spin current density $j_s$ in bulk iron at different depths, for different absorption profiles: Homogeneous absorption profile exciting only primary electrons (yellow curve), Lambert-Beer absorption profile with only primary excitations (black curve), and including primary and secondary electrons (red curve). The values for $j_s$ were obtained using eq. \eqref{eq::Spcrr}.}
    \label{fig:MaxJ_iron}

\end{figure}

As shown in fig. \ref{fig:MaxJ_iron}, contrary to other curves, where the maximum values increase until a certain depth, and then gradually decrease (red and black curves), the maximum value of the spin current density $j_s$ in a non-interactive picture (yellow curve), where only elastic scatterings are taking place, increases steadily in depth. This trend can be explained by examining the dynamics of particles following laser absorption. In the non-interactive picture (yellow curve), particles are initially excited based on a constant absorption profile throughout the material. However, as they move towards the depth of the material, the distribution of particles changes since their directionality is initially random. It is likely that particles will move ballistically in the absence of further interactions (change in direction or generation of new particles), and that the density of particles within the material will not be significant. Due to their initial distribution when they are created as well as their interactions with each other, particles, over time, acquire a common directionality through interaction. Hence, the absorption profile is responsible for the observed tendency in the maximum spin current density values in bulk iron at different depths.

\subsection{Spin current density in Fe/Au bilayer system}

We now proceed to examine the spin transport in a bilayer system containing iron and gold. Unlike bulk iron, laser radiation induces primary excitations within both layers in accordance with the absorption profile and these are initially distributed within the material accordingly, as illustrated in fig. \ref{fig:Prob_depth}. Fig. \ref{fig:Prob_depth} shows the data of the laser light absorption probability for different thicknesses of iron layers, which  determine the position of the interface in the bilayer system: $\SI{7}{\nano\meter}$ (yellow curve), $\SI{10}{\nano\meter}$ (blue curve) and $\SI{16}{\nano\meter}$ (black curve). These curves were obtained following the procedure described in ref.\cite{Eschenlohr,Khorsand2014} with absorption coefficients extracted from ref.\cite{Wolfgang2009}. The reason behind choosing these specific sizes is to study how smaller or larger iron layers
influence the dynamics. The data in fig. \ref{fig:Prob_depth} was calculated for a semi-infinite system, to avoid the reflective effects from gold.

\begin{figure}[H]
  \includegraphics[scale =0.6, center]{graphs2/Abs_intensity.eps}

  \caption{Absorption probability with respect to the depth for Fe/Au bilayer system with an interface at $\SI{7}{\nm}$ (yellow), $\SI{10}{\nm}$ (blue) and $\SI{16}{\nm}$ (black). Both materials have a different absorption profile. At the interface, one can observe a sudden increase in the absorption when the energy enters into the gold layer.}

\label{fig:Prob_depth}
\end{figure}

Over time, primary electrons generate secondary electrons, and these subsequently generate further electrons, a process that is repeated in both materials. At the boundary between iron and gold, which will be later addressed as interface, particles are transmitted in different proportion due to the spin dependent transmission probability (spin filtering effects).

In light of our previous analysis of iron, we now have a better understanding of how the absorption profile influences the behavior of the spin current density. Fig. \ref{fig:Max_J} illustrates how the absorption profile and thickness of iron layers affect the maximum value of spin current density in a bilayer system before and after crossing the interface. The magnitude of the peak of the spin current density is compared between the Lambert-Beer absorption profile in bulk iron (red curve), in the Fe/Au bilayer system (yellow curve) and using a calculated absorption profile (blue curve) (using parameters, procedure and absorption profile from refs. \cite{Wolfgang2009,Eschenlohr,Khorsand2014}, respectively).

\begin{figure}[H]
\centering
    \includegraphics[scale =0.7, center]{graphs2/Maxj_thickness_all.eps}

    \caption{A comparison of the maximum values of the spin current density $j_s$, obtained from eq. \eqref{eq::Spcrr}, for different thicknesses of iron. The value is compared between the Lambert-Beer absorption profile in bulk iron (red) and in the Fe/Au bilayer system (yellow), as well as one simulation using a calculated absorption profile (blue), applying eqs. from refs. \cite{Khorsand2014,Eschenlohr} and parameters in ref. \cite{Wolfgang2009}.}
    \label{fig:Max_J}

\end{figure}

In fig.~\ref{fig:Max_J} Iron (red) and Fe/Au (yellow) exhibit a large difference in maximum spin-current density as a result of the interface, which increases the value of $j_s$. However, the interface does not alter the tendency observed over a preferential thickness of iron layer. However, when the calculated absorption profile (from fig. \ref{fig:Prob_depth}) is taken, the preferred thickness changes and the larger spin current injection from iron to gold appears at a lower iron thickness. By analyzing the Lambert-Beer (yellow) and the calculated (blue) absorption profiles, one realizes the preferable thickness changes. This means that because of the influence of light absorption in the bilayer system, there is a specific thickness at which a larger spin current injection from iron to gold is observed. Spin dependent particles that cross between layers reconfigures the distribution of particles in gold, which influences the spin current density. This can be analyzed using the equation for spin polarization $P$ :
\begin{equation}\label{eq:spin_pol}
    P = \frac{\eta^{\uparrow} - \eta^{\downarrow}}{\eta^{\uparrow} + \eta^{\downarrow}},
\end{equation}
where $\eta^{\uparrow}$ ($\eta^{\downarrow}$) is the particle density of spin up (down) particles. Fig. \ref{fig:SpCrr_pol} shows the spin current density $j_s$ (top graph) and the spin polarization $P$ (lower graph), calculated using eqs. \eqref{eq::Spcrr} and \eqref{eq:spin_pol}, respectively, for three different thicknesses of iron layers.

\begin{figure}[H]
    \centering
    \includegraphics[scale =0.66, center]{graphs2/SpCrr_pol.eps}

    \caption{Spin current density $j_s$ (top graph) and spin polarization $P$ (bottom graph), calculated before (dashed lines) and after (solid lines) crossing the interface. The different thicknesses of the iron layers are displayed as color degraded lines.}

    \label{fig:SpCrr_pol}

\end{figure}

The dashed lines indicate quantities taken before crossing the interface, while the solid lines indicate quantities taken after crossing the interface. The lines depicting the different thicknesses of the iron layer are color degraded. For the spin current density, the observed behavior is similar to that of fig.~\ref{fig:etavel_J_BL}. In this bilayer system, we focus on the thickness of the iron layer instead of the propagation depth. As shown in fig.~\ref{fig:SpCrr_pol}, the maximum value of the spin current density $j_s$ increases up to a certain thickness (here $\SI{16}{\nano\meter}$), but decreases for a thicker film. As a result, it appears that the maximum value of spin current density does not increase with increasing iron layer thickness. Instead, there exists an optimal size of iron layer at which spin injection from iron to gold is more effective. The violet curves in fig.~\ref{fig:SpCrr_pol} are the spin polarization curves, which were calculated using eq. \eqref{eq:spin_pol}. They indicate that in time, always more particles with spin up will cross the interface as a result of the difference in transmission probability for the spins. By increasing the size of the iron layer, the maximum value of spin polarization will be reached at later times. This delay is a consequence of the time necessary for particles to interact within the material and travel toward the interface. As the thickness of iron layer increases, spin polarization increases since fewer particles cross over to gold when the layer is larger.

\section{Summary}

In conclusion, we investigated how different iron thicknesses affect the generation of spin currents at the interface of a Fe/Au bilayer system.  We first examined the effects of the light absorption profile on the spin current density within a bulk iron layer, and we found that the absorption profile influences the spin current density non monotonously.
For the Fe/Au bilayer, we computed the nonequilibrium transport of excited electrons and studied how different absorption profiles influence the spin current injection from one layer into another.
This analysis also revealed that secondary electron generation and the interface affect spin transport only in terms of magnitude. We then investigated the combined effect of light absorption and spin-dependent interface transmission probability  on the spin current density in the Fe/Au bilayer system.
The magnitude of the spin current density from Fe to Au is influenced by the thickness of the Fe layer and shows a maximum at intermediate thicknesses. We found that the absorption characteristics of the exciting laser pulse has a direct influence on the
efficiciency of the spin injection from one layer into another.

\ack Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - TRR 173/2 - 268565370 Spin+X (Project B03). B.S. acknowledges financial support by the Dynamics and Topology Center, funded by the State of Rhineland-Palatinate.

\section*{References}

\bibliographystyle{iopart-num}