\section{Introduction}
\IEEEPARstart{I}{n} Part I of this two parts paper, we presented the advantages, the reference architecture and three core applications of wireless federated learning (WFL), the main challenges toward its efficient integration in the sixth generation of wireless networks (6G), and the identified future directions \cite{part1}. This second  part is motivated by the fact that, as it has  been thoroughly discussed in II. A and III.A of the first part, the utilization of advanced multiple access protocols and the joint optimization of the communication and computing resources can facilitate WFL to meet the stringent latency requirements of 6G networks \cite{part1, konevcny2016, 6G}. To this direction, the use of non-orthogonal multiple access (NOMA) has been proposed in the first part, due to offering low-latency and improved fairness by serving multiple users in the same resource block, in opposition to  orthogonal multiple access (OMA), where each user occupies a single resource block \cite{ding2017}.
\par  The aim of the second part of this paper is to investigate the potential of using NOMA in WFL, by focusing on a specific paradigm. In more detail, we introduce the \textit{Compute-then-Transmit NOMA (CT-NOMA)} protocol, which is based on the use of two phases during a WFL round. Specifically, users firstly execute and complete the local model training concurrently and  afterwards upload the trained model to the server simultaneously via NOMA. Moreover, we consider two detection schemes for the implementation of  successive interference cancellation (SIC) process in NOMA that correspond to the utilization of a fixed or variable decoding order. The later, is implemented through the time-sharing (TS) strategy, which improves the performance of uplink NOMA by achieving any point of the multiple access channel (MAC) capacity region \cite{goldsmith}. Furthermore, we minimize the total delay of a communication round subject to users' energy constraints for both considered schemes, by jointly optimizing the available computation and communication resources, i.e., users' CPU frequency for local training, the transmission energy and the time intervals for local computations and parameter transmission.  Finally, simulation results demonstrate the effectiveness of the proposed protocol in reducing the delay of a WFL communication round, compared the considered benchmark that is based on time-division multiple access (TDMA).
\vspace{-0.1in}
\section{System Model and Proposed Protocols}
\subsection{Wireless federated learning model}
Based on the reference architecture that has been presented in detail in the first part of this work, we consider a WFL learning system which consists of $N$ users indexed as $n \in \mathcal{N}=\{1,2,...,N\}$ and a server/base station (BS). Each user $n$ has a local dataset $\mathcal{D}_n$,  where $D_n=\vert \mathcal{D}_n \vert$ are the total data samples. Next we briefly describe the steps of an arbitrary, $i$-th communication round, between the users and the BS in the WFL system.
\begin{enumerate}[i)]
	\item The BS broadcasts the global parameter $\boldsymbol{w}^i$ to all users during the considered round.
	\item After receiving the global model parameter, each user $n \in \mathcal{N}$, trains the local model through its dataset, and then uploads the trained local parameter $\boldsymbol{w}^{i+1}_n$ to the server.
	\item After receiving all the local parameters, the server aggregates them, in order to update the global model parameter  $\boldsymbol{w}^{i+1}$.
\end{enumerate}
The whole procedure is repeated until the global model converges. Moreover, during the first round the server initializes  $\boldsymbol{w}^0$.
\subsection{Compute-then-Transmit NOMA}
As it has already been mentioned, CT-NOMA is based on the use of two consecutive phases during a WFL round, i.e., the computation and the communication phase. Specifically, in the first phase users  execute the local computations, while in the second phase they transmit their messages, i.e., the trained parameters, to the BS. During the information transmission phase, the NOMA protocol is considered. According to NOMA, users are capable of transmitting their messages simultaneously, while the BS decodes the users' messages by utilizing SIC. As implied by the term CT-NOMA, all users begin the information transmission procedure via NOMA at the same time instant. Therefore, all users are constrained to  complete the local computations before the information transmission phase, i.e., within time duration  $\tau$, while the transmission phase duration is denoted as $t$. The CT-NOMA protocol is illustrated in Fig. \ref{CTNOMA}.
\par
The utilized computation resources for local model training,
i.e., CPU cycle frequency, from the $n$-th user is denoted as $f_n$. The number of CPU cycles for the $n$-th user to perform
one sample of data in local model training is denoted by $c_n$. Hence, the computation time dedicated for a local iteration is given as
\begin{equation}
\tau_n=\frac{c_nD_n}{f_n}, \quad \forall n \in \mathcal{N}.
\end{equation}
 Accordingly, the energy consumption for a local iteration, can be expressed as follows
\begin{equation}
E^{\mathrm{comp}}_n=\zeta c_n D_nf^2_n = \zeta \frac{c^3_nD^3_n}{\tau_n^2}, \quad \forall n \in \mathcal{N},
\end{equation}
where $\zeta$  is a constant parameter related to the hardware
architecture of device $n$.
As it was discussed previously, all users are enforced to complete the local computations within $\tau$, with the corresponding energy that is consumed by each user being a decreasing function with respect to $\tau$.
Thus, it should hold
\begin{equation}
	\tau_n=\tau, \quad \forall n \in \mathcal{N}.
\end{equation}

\par
 In the continue, the two proposed schemes regarding information transmission will be described, termed as NOMA with fixed decoding order (FDO)  and NOMA with TS.
\subsubsection{NOMA with FDO} Fixed decoding order refers to the classical uplink NOMA protocol with SIC, which is widely used in the literature \cite{ding2017}. Let $\pi_n$ be the the $n$-th user's decoding order position at the BS, while without loss of generality we consider that $\pi_1<\pi_2<...<\pi_N$. According to NOMA, for decoding the first user's message, interference is created by the residual users, i.e., $n=2,...,N$, while on the second user's message, interference is created by the users indexed as $n=3,...,N$ and so on. Finally the last user to be decoded, user $N$, experiences no interference. Following that, the data size $Z_n$ of the $\boldsymbol{w}_n$ model parameter, that the $n$-th user transmits within $t$ time duration, should satisfy the following condition
\begin{equation}
Z_n \leq tB\log_2\left(1 + \frac{\frac{E_n}{t}g_n}{\sum_{\{i \in \mathcal{N}|\pi_i>\pi_n\}}\frac{E_i}{t}g_i+BN_0}\right),  \forall n \in \mathcal{N},
\end{equation}
where $B$ is the available bandwidth, $N_0$ denotes the power spectral density, while $E_n$ is the $n$-th user's consumed energy for the aforementioned transmission, with $\frac{E_n}{t}$ being the corresponding  transmit power. Moreover, $g_n=\vert h_n \vert^2d_n^{-2}$ denotes the channel gain, where the complex random variable $h_n \sim \mathcal{CN}(0,1)$
is the small scale fading and $d_n$ is the distance between user
$n$ and the BS. Finally, we assume that the channel remains constant during a communication round and can be perfectly estimated by the BS.
\subsubsection{NOMA with TS} The basic principle of this scheme is that the users' decoding order can change throughout the $t$ transmission time duration, unlike to the NOMA with FDO where a fixed decoding order is considered. According to \cite{goldsmith}, the capacity region which can be achieved by uplink NOMA with TS, known as multiple access channel (MAC) capacity region, is described as \cite{goldsmith}
 \begin{equation}
\sum_{i \in \mathcal{S}}^{}Z_{i} \leq tB\log_2\left(1 + \frac{\sum_{i \in \mathcal{S}}^{}E_{i}g_i}{tBN_0}\right), \quad \forall \mathcal{S} \subseteq \mathcal{N}, \, \mathcal{S} \neq {\O},
\label{MACcr}
\end{equation}
where $\mathcal{S}$ is a non-empty subset of $\mathcal{N}$. Furthermore, it deserves to be mentioned that \eqref{MACcr} also determines the performance of  rate splitting multiple access (RSMA), which is an alternative scheme to achieve any point of the MAC capacity region \cite{rimoldi1996,RSMA}.
\begin{figure}[h!]
\vspace{-0.1in}
	\centering
	\includegraphics[width=0.8\linewidth]{protocol.pdf}
	\vspace{-0.1in}
	\caption{\textit{Compute-then-Transmit NOMA} during a WFL round.}
	\label{protocol}
	\vspace{-0.2in}
	\label{CTNOMA}
\end{figure}
\section{Delay Minimization  of a WFL Round}
The aim of the considered optimization problem is to minimize the total delay of a WFL round, i.e., $T=\tau + t$, which is the sum of the computation and transmission latency.  Note that, since the transmit power of the BS is much higher than that of
the users’ and the BS transmits the same message to all users, we ignore the delay of the server for broadcasting
the global parameter. Hereinafter, we focus on a single arbitrary round, while the subsequent analysis can be similarly carried out for each round. Finally, we assume that each user transmits the same data size $Z_n=Z, \, \forall n \in \mathcal{N}$, related with the model parameters, while the considered data size can be predetermined by the central server.
\par
The maximum available CPU clock speed of each user is denoted as $f^{\mathrm{max}}_n$. Thus, it should hold $f_n \leq f^{\mathrm{max}}_n, \,  \forall n \in \mathcal{N}$, which is equivalent to $\tau \geq \underset{n \in \mathcal{N}}{\mathrm{max}}\left(\frac{c_nD_n}{f^{\mathrm{max}}_n}\right)\triangleq a_1$. Moreover, the maximum available energy of each user is $E^{\mathrm{max}}_n$. Hence, it should hold
\begin{equation}
	E^{\mathrm{comp}}_n + E_n=\zeta \frac{c^3_nD^3_n}{\tau^2}+E_n \leq E^{\mathrm{max}}_n, \quad \forall n \in \mathcal{N},
\end{equation}
since the total consumed energy  for both computation and communication purposes, cannot exceed the maximum available energy.
\subsection{CT-NOMA with time-sharing}
The optimization problem for minimizing the total latency in the case of CT-NOMA with TS, can be written as
\begin{equation}
\begin{aligned}
& &  &\underset{\text{\emph{$\tau,t$},{\boldmath{$E$}}}}{\textbf{min}} \quad
\text{\emph{{$\tau+t$}}} \\
& & &\textbf{\emph{s.t.}}
\quad \mathrm{C}_{1}: 	\zeta \frac{c^3_nD^3_n}{\tau^2}+E_n \leq E^{\mathrm{max}}_n, \quad \forall n \in \mathcal{N}, \\
& & &\qquad \mathrm{C}_{2}: 	\sum_{i \in \mathcal{S}}^{}Z \leq tB\log_2\left(1 + \frac{\sum_{i \in \mathcal{S}}^{}E_ig_i}{tBN_0}\right), \quad \forall \mathcal{S} \subseteq \mathcal{N}, \\
& & & \qquad   \mathrm{C}_3: \tau \geq a_1, \,t\geq 0,E_n \geq 0, \, \forall n \in \mathcal{N},
\end{aligned}\label{primary}
\end{equation}
where  $\mathrm{C}_{2}$ is related to the users' achievable capacity region in \eqref{MACcr}. It is easy to verify that the problem in \eqref{primary} is jointly convex with respect to (w.r.t.) $\tau,t,\boldsymbol{E}$ and thus can optimally be solved with standard convex-optimization solving methods, such as the interior-point. However, the computational cost to directly solve this by standard methods remains intractable, due to the large number of constraints in $\mathrm{C}_{2}$. More specifically, the number of constraints in $\mathrm{C}_{2}$ is equal to $2^N-1$, since this is the total number of all non-empty subsets $\mathcal{S}$ of the set $\mathcal{N}$. As a result, solving \eqref{primary} for large $N$ comes at the expense of high complexity. To alleviate this burden, an efficient method for solving \eqref{primary} will be proposed. Firstly, we consider a fixed value of $\tau$. Since $\log(\cdot)$ is an ascending  function w.r.t $E_n$, an optimal solution occurs when the constraints $\mathrm{C}_{1}$ are satisfied with equality, i.e., users utilize their whole available energy, which leads to
\begin{equation}\label{energy_eq}
E_n=E^{\mathrm{max}}_n-\zeta \frac{c^3_nD^3_n}{\tau^2}, \quad  \forall n \in \mathcal{N}.
\end{equation}
Furthermore, since $E_n \geq 0, \,  \forall n \in \mathcal{N}$, by manipulating \eqref{energy_eq}, yields $\tau \geq \underset{n \in \mathcal{N}}{\mathrm{max}}\left(\sqrt{\zeta\frac{c^3_nD^3_n}{E^{\mathrm{max}}_n}}\right) \triangleq a_2$. Thus, by also recalling that $\tau \geq a_1$, it should finally hold for $\tau$
\begin{equation}
\tau \geq \mathrm{max}\{a_1,a_2\}\triangleq \tau_{\mathrm{low}}.
\end{equation}
Next, after substituting \eqref{energy_eq} and considering a known $\tau \geq \tau_{\mathrm{low}}$, the optimization problem in \eqref{primary} can be reformulated as
\begin{equation}
\begin{aligned}
& &  &\underset{\text{\emph{$t$}}}{\textbf{min}} \quad
\text{\emph{{$t$}}} \\
& & &\textbf{\emph{s.t.}} \quad	\sum_{i \in \mathcal{S}}^{}Z \leq tB\log_2\left(1 + \frac{\sum_{i \in \mathcal{S}}^{}A_i}{tBN_0}\right), \quad \forall \mathcal{S} \subseteq \mathcal{N}, \\
\end{aligned}\label{t_primary_problem}
\end{equation}
where $A_i$ is given by
\begin{equation}
	A_i=E_ig_i=\left(E^{\mathrm{max}}_i-\zeta \frac{c^3_iD^3_i}{\tau^2}\right)g_i.
\end{equation}
Note that $A_i$, is constant for a given $\tau$.
Following that, without loss of generality we assume that $A_1\geq A_2\geq...\geq A_N$. According to \cite{diamantoulakis}, when NOMA is combined with the TS scheme and by considering the aforementioned order of $A_i, \, \forall i \in \mathcal{N}$, the achievable capacity region of the users can be equivalently described as
\begin{equation}
\sum_{i=n}^{N}Z \leq tB\log_2\left(1 + \frac{\sum_{i=n}^{N}A_i}{tBN_0}\right), \quad \forall n \in \mathcal{N},
\label{cr1}
\end{equation}
while by considering that $\sum_{i=n}^{N}Z=Z(N+1-n)$, \eqref{cr1} can be rewritten as
\begin{equation}
	Z \leq \frac{tB\log_2\left(1 + \frac{\sum_{i=n}^{N}A_i}{tBN_0}\right)}{(N+1-n)}, \quad \forall n \in \mathcal{N}.
	\label{cr2}
\end{equation}
It should be noted that the capacity region is now described by $N$ constraints, i.e., $Z$ is bounded by a set of $N$ inequalities.
Next, by exploiting the alternative, but still equivalent representation of the capacity region in \eqref{cr2}, the optimization problem in \eqref{t_primary_problem} can be rewritten as
\begin{equation}
\begin{aligned}
& &  &\underset{\text{\emph{$t$}}}{\textbf{min}} \quad
\text{\emph{{$t$}}} \\
& & & \textbf{\emph{s.t.}} \quad	Z(N+1-n) \leq tB\log_2\left(1 + \frac{\sum_{i=n}^{N}A_i}{tBN_0}\right),
 \forall n \in \mathcal{N}.\\
\end{aligned}\label{t_problem}
\end{equation}
\begin{lemma}
	The optimal value of $t$ for the problem in \eqref{t_problem} can be written as
	\begin{equation}\label{t_solution}
	t^*=\underset{n \in \mathcal{N}}{\mathrm{max}}\left[-\frac{z_n\ln(2)\sum_{i=n}^{N}A_i}{B\left(z_nN_0\ln(2)+\mathcal{W}_{-1}(b_n)\sum_{i=n}^{N}A_i\right)}\right],
	\end{equation}
	where $z_n=(N+1-n)Z$ and $b_n$ is given by
	\begin{equation}
	b_n=-\frac{2^{-\frac{z_nN_0}{\sum_{i=n}^{N}A_i}}z_nN_0\ln(2)}{\sum_{i=n}^{N}A_i},
	\end{equation}
	while $\mathcal{W}_{-1}(\cdot)$ denotes the secondary branch of the Lambert W function.
\end{lemma}
\begin{IEEEproof}
	 We will show that the optimal $t$ is given by the most stringent inequality, among the set of $N$ inequalities in \eqref{t_problem}, when this hold with equality. Firstly, it is straightforward to show that the function
	\begin{equation}
		f(t)=tB\log_2\left(1 + \frac{\sum_{i=n}^{N}A_i}{tBN_0}\right), \quad \forall t>0,
	\end{equation}
	is monotonically increasing with respect to $t$. Following that, we assume that $t'$ is optimal and satisfies
	\begin{equation}
		Z(N+1-n) < t'B\log_2\left(1 + \frac{\sum_{i=n}^{N}A_i}{t'BN_0}\right), \, \forall n \in \mathcal{N},
	\end{equation}
	i.e., all constraints are satisfied with strict inequality.
	It will be shown by contradiction that this is not an optimal solution for $t$.
	Since $f(t)$ is increasing w.r.t $t$, there exist an $t^*<t'$, for which at least one inequality constraint from \eqref{t_problem} is satisfied with equality. This observation contradicts the fact that $t'$ is optimal, since the objective is to minimize $t$. Therefore, $t^*$ is given by the most stringent inequality, among the set of $N$ inequalities in \eqref{t_problem}, when this hold with equality.
	After some mathematical manipulations, $t^*$ can be written as in \eqref{t_solution} and the proof is completed.
\end{IEEEproof}
As a matter of fact, the problem can be efficiently solved with closed form solutions for a fixed $\tau$. In the continue, we will show that the global optimal solution can be obtained via the bisection method, in order to find the optimal value of $\tau$. Firstly, it will be shown that $\tau^*$ is bounded, according to the following lemma.
\begin{lemma}
	The optimal $\tau$ is bounded in the following interval
	\begin{equation}
		\tau_{\mathrm{low}} \leq \tau^* \leq \tau_{\mathrm{up}},
	\end{equation}
where $\tau_{\mathrm{up}}=T(\bar{\tau})=\bar{\tau}+t^*(\bar{\tau})$ is the total delay of the communication round that corresponds to an arbitrary feasible value of $\tau$, denoted by $\bar{\tau}$, while  $t^*(\bar{\tau})$ denotes the optimal $t$ for a fixed $\tau = \bar{\tau}$ and is given by \eqref{t_solution}.
\end{lemma}
\begin{IEEEproof}
 As mentioned previously, the optimization problem is feasible when the condition $\tau \geq \tau_{\mathrm{low}}$ is satisfied.  Since $\bar{\tau}$ could be any feasible solution but not necessarily the optimal one, it holds $T(\bar{\tau}) \geq T(\tau^*)$. Following this, by considering that the dedicated time for computations cannot exceed the total delay for any feasible $\tau$, i.e., $\tau < T(\tau)$, we conclude to $\tau^* < T(\tau^*) \leq T(\bar{\tau})$. Thus, the upper bound of $\tau^*$ can be expressed as $\tau_{\mathrm{up}}=T(\bar{\tau})$ and the proof is completed.
\end{IEEEproof}
\par
 As a result, we have derived the upper and lower bound of $\tau^*$ and now we are ready to apply the bisection method in the considered interval and propose Algorithm 1 for obtaining the optimal solutions that minimize the delay of a communication round. Next, the main steps of the proposed algorithm are discussed. After the necessary initializations in line 1, the bisection method is applied throughout lines 2-13. More specifically, in lines 3- 6, the total delay of the communication round is calculated at the points $\tau=\tau_{\mathrm{m}}$ and $\tau=\tau_{\mathrm{up}}$.  Following that, in lines 7-12, by comparing $T(\tau_{\mathrm{m}})$ and $T(\tau_{\mathrm{up}})$, the bounds of $\tau$ are properly adjusted in each iteration, until the convergence is achieved. It should be highlighted that the joint convexity of the primary problem in \eqref{primary}, w.r.t. all the considered variables, guarantees the convergence of the bisection method to the optimal solution. After the resolution of the optimization problem and the calculation of $\tau^*,t^*$, the optimal $E_n,f_n$ can be given as $E^*_n=E^{\mathrm{max}}_n-\zeta \frac{c^3_nD^3_n}{\tau^{*2}}$ and $f^*_n=\frac{c_nD_n}{\tau^*}$. To this end, the major complexity of Algorithm 1 lies in applying the bisection method and in searching among $N$ values, in order to derive $t^*$ from \eqref{t_solution}. As a result, the complexity can be expressed as $\mathcal{O}\left(N\log_2\left(\frac{\tau_{\mathrm{up}}-\tau_{\mathrm{low}}}{\epsilon}\right)\right)$, where $\epsilon$ denotes the algorithm's tolerance.
\makeatletter
\def\State\hskip-\ALG@thistlm{\State\hskip-\ALG@thistlm}
\makeatother
\begin{algorithm}
	\caption{Delay Minimization for CT-NOMA with TS}
	\begin{algorithmic}[1]

		\State\hskip-\ALG@thistlm  \textbf{Initialize} \text{$\tau_{\mathrm{low}}=\mathrm{max}\{a_1,a_2\},\,\tau_\mathrm{m}=\bar{\tau},\,\tau_{\mathrm{up}}=T(\bar{\tau}),\,\epsilon;$}
		\State\hskip-\ALG@thistlm \textbf{while} \text{$\tau_{\mathrm{up}}-\tau_{\mathrm{low}}>\epsilon$} \textbf{do}
		\State\hskip-\ALG@thistlm \quad \textbf{Set} $\tau=\tau_\mathrm{m},$ \text{ and derive} $t^*(\tau_\mathrm{m})$ \text{from \eqref{t_solution}}
		\State\hskip-\ALG@thistlm \quad \textbf{Set} $T(\tau_\mathrm{m})=\tau_\mathrm{m}+t^*(\tau_\mathrm{m});$
		\State\hskip-\ALG@thistlm \quad \textbf{Set} $\tau=\tau_\mathrm{up},$ \text{and derive} $t^*(\tau_\mathrm{up})$ \text{from \eqref{t_solution}}
		\State\hskip-\ALG@thistlm \quad \textbf{Set} $T(\tau_\mathrm{up})=\tau_\mathrm{up}+t^*(\tau_\mathrm{up});$
		\State\hskip-\ALG@thistlm \quad \textbf{if} $T(\tau_\mathrm{m})<T(\tau_\mathrm{up})$
		\State\hskip-\ALG@thistlm \qquad $\tau_\mathrm{up}=\tau_\mathrm{m};$
		\State\hskip-\ALG@thistlm \quad \textbf{else}
		\State\hskip-\ALG@thistlm \qquad $\tau_\mathrm{low}=\tau_\mathrm{m};$
		\State\hskip-\ALG@thistlm \quad \textbf{end if}
		\State\hskip-\ALG@thistlm \quad $\tau_\mathrm{m}=\frac{\tau_\mathrm{low}+\tau_\mathrm{up}}{2};$
		\State\hskip-\ALG@thistlm \textbf{end while}
		\State\hskip-\ALG@thistlm \textbf{Output} $\tau^* =\tau_\mathrm{up}, \, t^*=t^*(\tau_\mathrm{up}),\, T^*=\tau^*+t^*;$
	\end{algorithmic}\label{alg}
\end{algorithm}

\subsection{CT-NOMA with fixed decoding order}
The optimization problem for minimizing the total delay in the case of CT-NOMA with FDO, can be written as
\begin{equation}
\begin{aligned}
& &  &\underset{\text{\emph{$\tau,t$},{\boldmath{$E,\pi$}}}}{\textbf{min}} \quad
\text{\emph{{$\tau+t$}}} \\
& & &\textbf{\emph{s.t.}}
\quad \mathrm{C}_{1}: 	\zeta \frac{c^3_nD^3_n}{\tau^2}+E_n \leq E^{\mathrm{max}}_n, \quad \forall n \in \mathcal{N}, \\
& & &\qquad \mathrm{C}_{2}: 	Z \leq tB\log_2\left(1 + \frac{\frac{E_n}{t}g_n}{\sum_{\{i \in \mathcal{N}|\pi_i>\pi_n\}}\frac{E_i}{t}g_i+BN_0}\right)\\
 & & &\qquad \qquad \forall n \in \mathcal{N}, \\
& & & \qquad   \mathrm{C}_3: \tau \geq a_1, \,t\geq 0,E_n \geq 0, \, \forall n \in \mathcal{N}, \, \boldsymbol{\pi}\in \Pi,
\end{aligned}\label{fixed_decoding}
\end{equation}
where the permutation $\boldsymbol{\pi}$ belongs to the set $\Pi$, defined as the set of all possible decoding orders among the $N$ users.
 Note that the total number of permutations, i.e., all possible decoding orders, is equal to $\vert\Pi\vert=N!$. Therefore, the optimization problem in \eqref{fixed_decoding} is of a combinatorial nature. However, an exhaustive search among all possible decoding orders is prohibitive, due to high complexity. To this end, we adopt an ascending decoding order w.r.t. to the channel gains, according to which the messages of the users with  weaker channel gains are exposed to less interference during the decoding process. It is notable that this a common selection for uplink NOMA systems, since it provides fairness among users without reducing the sum rate \cite{yang2016,diamantoulakis}. Accordingly, by considering  without loss of generality that $g_1\geq g_2\geq...\geq g_N$, we set $\pi_1<\pi_2<...<\pi_N$. As a result, $\mathrm{C}_{2}$ of \eqref{fixed_decoding} can be reformulated as
\begin{equation}
	Z \leq tB\log_2\left(1 + \frac{E_ng_n}{\sum_{i=n+1}^{N}E_ig_i+tBN_0}\right), \quad \forall n \in \mathcal{N}.
	\label{C2initial}
\end{equation}
By substituting \eqref{C2initial} in \eqref{fixed_decoding}, the later can be written as
\begin{equation}
\begin{aligned}
& &  &\underset{\text{\emph{$\tau,t$},{\boldmath{$E$}}}}{\textbf{min}} \quad
\text{\emph{{$\tau+t$}}} \\
& & &\textbf{\emph{s.t.}}
\quad \mathrm{C}_{1}: 	\zeta \frac{c^3_nD^3_n}{\tau^2}+E_n \leq E^{\mathrm{max}}_n, \quad \forall n \in \mathcal{N}, \\
& & &\qquad \mathrm{C}_{2}: Z \leq tB\log_2\left(1 + \frac{E_ng_n}{\sum_{i=n+1}^{N}E_ig_i+tBN_0}\right),\\
& & &\qquad \qquad \forall n \in \mathcal{N}, \\
& & & \qquad   \mathrm{C}_3: \tau \geq a_1, \,t\geq 0,E_n \geq 0, \, \forall n \in \mathcal{N},
\end{aligned}\label{fixed_decoding2}
\end{equation}
which is non-convex due to the coupling of $t$ and $\boldsymbol{E}$ in constraint $\mathrm{C}_{2}$. In order to address the non-convexity issues, we first rewrite $\mathrm{C}_{2}$ as
\begin{equation}
	2^{\frac{Z}{tB}} \leq \frac{E_ng_n}{\sum_{i=n+1}^{N}E_ig_i+tBN_0}, \, \forall n \in \mathcal{N}.
\end{equation}
In the continue we set $E_n\triangleq \exp(\tilde{E}_n),\, \forall n \in \mathcal{N}$ and $t\triangleq \exp(\tilde{t})$. After some mathematical manipulations, $\mathrm{C}_{2}$ can be written as
\begin{equation}
\begin{aligned}
& & &	\ln\left(2^{\frac{Z}{B}\exp(-\tilde{t})}-1\right)\\
& & &+\ln\left(\sum_{i=n+1}^{N}\exp(\tilde{E}_i-\tilde{E}_n)g_i+\exp(\tilde{t}-\tilde{E}_n)BN_0\right)\\
& & & \leq \ln(g_n), \quad \forall n \in \mathcal{N}.
\end{aligned}
\end{equation}
By exploiting the aforementioned formulation of $\mathrm{C}_{2}$, the optimization problem in \eqref{fixed_decoding2} can be transformed to the following equivalent one:
\begin{equation}
\begin{aligned}
& &  &\underset{\text{\emph{$\tau,\tilde{t}$},{\boldmath{$\tilde{E}$}}}}{\textbf{min}} \quad
\text{\emph{{$\tau+\exp(\tilde{t})$}}} \\
& & &\textbf{\emph{s.t}}
\quad \mathrm{C}_{1}: 	\zeta \frac{c^3_nD^3_n}{\tau^2}+\exp(\tilde{E}_n) \leq E^{\mathrm{max}}_n, \quad \forall n \in \mathcal{N}, \\
& & &\qquad \mathrm{C}_{2}: 	\ln\left(2^{\frac{Z}{B}\exp(-\tilde{t})}-1\right)
+\ln\Biggl(\sum_{i=n+1}^{N}\exp(\tilde{E}_i-\tilde{E}_n)g_i\\
& & & \qquad \qquad +\exp(\tilde{t}-\tilde{E}_n)BN_0\Biggr) -\ln(g_n) \leq 0, \, \forall n \in \mathcal{N}, \\
& & & \qquad   \mathrm{C}_3: \tau \geq a_1, \,\tilde{t}\in \mathbb{R},{\tilde{E}_n} \in \mathbb{R}.
\end{aligned}\label{fixed_decoding3}
\end{equation}
It can be easily shown that the optimization problem in \eqref{fixed_decoding3} is convex, while the proof is omitted for brevity. Thus, it can efficiently be solved with standard convex optimization methods.
\section{Performance Evaluation and Discussion}
Table I summarizes the simulations' parameters setup. The following figures have been extracted by means of Monte Carlo.
\begin{figure}[h!]
	\centering
	\includegraphics[width=0.8\linewidth]{Latency_vs_energy.pdf}
	\vspace{-0.1in}
	\caption{Impact of the users' maximum available energy on latency, with $Z=0.8\mathrm{Mbits}$.}
	\vspace{-0.1in}
	\label{latency}
\end{figure}
In Fig. \ref{latency}, the impact of the users' maximum available energy on the average latency during a WFL round is depicted. We use the TDMA-based protocol as benchmark  and compare it with the CT-NOMA protocol in terms of delay reduction, taking into account both considered detections schemes for NOMA. It can been observed that CT-NOMA with FDO slightly outperforms TDMA. However, CT-NOMA with TS clearly dominates both TDMA and CT-NOMA with FDO. This is due to the fact that CT-NOMA with TS can achieve any point of the capacity region, in contrast to CT-NOMA with FDO which achieves only the corner points. As a result, the efficient exploitation of the capacity region, as well as subsequently the flexible interplay among users' data rates that CT-NOMA with TS enables, lead to decreased delay during the WFL round. Moreover, in Fig. \ref{datasize}, the impact of the users' parameter data size on latency is illustrated. The superiority of CT-NOMA with TS against CT-NOMA with FDO, and TDMA is again corroborated. Furthermore, CT-NOMA with FDO presents a slightly enhanced performance in comparison with TDMA.
\begin{figure}[h!]
	\centering
	\includegraphics[width=0.8\linewidth]{Latency_vs_datasize.pdf}
		\vspace{-0.1in}
	\caption{Impact of the users' parameter data size on latency, with $E^{\mathrm{max}}_n=2\, \mathrm{Joule}$.}
		\vspace{-0.1in}
	\label{datasize}
\end{figure}
\par
By evaluating the demonstrated results, it is evident that CT-NOMA has the potential to provide decreased latency during a WFL round and subsequently accelerate the training process, which is an important requirement for the efficient integration of WFL in 6G. Moreover, since NOMA can achieve any point of the capacity region through the use of time-sharing (or rate splitting) technique, it can serve in future research as a performance upper bound in order to evaluate other multiple access schemes for WFL. Furthermore, the proposed protocol could serve as a baseline for more sophisticated network implementations, such as hybrid NOMA/OMA configurations, which could possibly enhance the scalability of WFL.
\begin{table}
	\caption{Simulation Settings}
	\centering
	\begin{tabular}{ p{1.5cm}|p{1.5cm}||p{1.5cm}|p{2cm}}
		\hline
		\centering
		\textbf{Parameter} & \textbf{Value} & \textbf{Parameter} & \textbf{Value}\\
		\hline\hline
		$f^{\mathrm{max}}_n$ & $1.5\mathrm{GHz}$ & $D_n$ & 1Mbit \\
		\hline
		$B$ & 1MHz & $N_0$ & 174dBm/Hz \\
		\hline
		$\zeta$ & $10^{-27}$ & $c_n$(cycles) & $\sim \mathcal{U}(10,40)$ \\
		\hline
		$N$ & 10 users & $d_n$ & $\sim \mathcal{U}(0,500\mathrm{m})$ \\
		\hline

	\end{tabular}
	\vspace{-0.2in}
\end{table}
	\vspace{-0.2in}

\bibliographystyle{IEEEtran}


\section{Introduction}
Ultrafast manipulation of electrons' spin remains at the heart of future generation spin-based memory technology \cite{bigot09,Stanciu2007,Kimel_2007_JPCM_review}. It has been observed that a fs laser pulse is capable of demagnetizing a ferromagnetic material~\cite{Bigot1996,Koopmans2000,Koopmans2003JPCM}. On the other hand, using these ultrashort pulses, magnetic switching has been reported in ferrimagnetic~\cite{Mangin2014,Hassdenteyfel2013,Kimel2005} and ferromagnetic materials~\cite{lambert14,John2017}. These observations have been explained through the spin dynamics within Landau-Lifshitz-Gilbert (LLG) equation of motion~\cite{Ostler2012,Wienholdt2013,Gerlach2017,Frietsch2020}.

The phenomenological LLG spin dynamics consists of spin precession and a transverse damping~\cite{landau35,Gilbert1955,Gilbert2004}. Such an equation of motion has been derived from a relativistic Dirac theory, where the transverse damping is found to originate from spin-orbit coupling~\cite{hickey09,Mondal2015a,Mondal2016,Mondal2018PRB}. However, at ultrashort timescales, the traditional LLG equation needs to be supplemented by several other spin torque terms~\cite{Mondal2019PRB}. Especially, at the ultrafast timescales, the magnetic inertia becomes particularly relevant~\cite{Henk2012}. The effect of magnetic inertia has been incorporated within extended LLG dynamics as a torque due to the second-order time derivative of the magnetization $\bm{M}(\bm{r},t)$. The inertial LLG (ILLG) equation of motion reads~\cite{Ciornei2010thesis,Ciornei2011,Bhattacharjee2012}
\begin{align}
    \frac{\partial \bm{M}}{\partial t} & = \bm{M} \times \left[-\gamma \bm{H} +  \frac{\alpha}{M_0}\frac{\partial \bm{M}}{\partial t} + \frac{\eta}{M_0} \frac{\partial^2 \bm{M}}{\partial t^2}\right]\,,
    \label{ILLG-eq}
\end{align}
where $M_0$ and $\bm{H}$ define the ground state magnetization and an effective field, respectively. The first and second terms in Eq.~(\ref{ILLG-eq}) represent the traditional LLG equation~\cite{Gilbert2004}.
The inertial spin dynamics in the last term of Eq.~(\ref{ILLG-eq}) gives rise to the spin nutation~\cite{Wegrowe2012,Wegrowe2016JPCM}. { The ILLG equation signifies the fact that the dynamics of a magnetic moment shows precession with nutation at ultrafast timescales, followed by transverse damping \cite{Henk2012}. The ILLG equation has schematically been depicted in Fig.~\ref{Fig1}. }
{ A simple dimension analysis shows that }the transverse damping is characterized by a dimensionless parameter $\alpha$, and the inertial dynamics are strengthened by inertial relaxation time $\eta$.
\begin{figure}[hbt!]
    \centering
    \includegraphics[scale = 1]{Fig1.pdf}
    \caption{ Schematic depiction of ILLG equation of motion. }
    \label{Fig1}
\end{figure}
The ILLG dynamics have been derived within the relativistic Dirac framework as well, where it shows that the Gilbert damping $\alpha$ and inertial relaxation time $\eta$ are tensors~\cite{Mondal2017Nutation}. In particular, the relativistic theory derives that the Gilbert damping dynamics is associated with the imaginary part of the susceptibility, while the inertial dynamics is given by the real part~\cite{Mondal2018JPCM}. Such findings are found to be consistent with a linear response theory of ferromagnet~\cite{Thonig_2014}. The inertial dynamics have also been derived within classical mechanics of a current loop~\cite{Giordano2020}. Eq.~(\ref{ILLG-eq}) has been applied to a single sublattice ferromagnet beyond ferromagnetic resonance (FMR), observing an additional peak due to nutation resonance~\cite{Olive2012,Olive2015,Makhfudz2020}. While the FMR peak appears at the GHz regime, the nutation resonance peak appears at the THz regime~\cite{cherkasskii2020nutation}. The ILLG equation has also been applied to antiferromagnets and ferrimagnets, and it has been predicted that the spin nutation should be better detected in antiferromagnets as it is exchange enhanced~\cite{Mondal2020nutation}.

Recently, the spin nutation resonance has been observed for ferromagnets in the experiment ~\cite{neeraj2019experimental}. Indeed, the nutation resonance peak has been seen at around 0.5 THz. Note that the experiment was performed in two-sublattice ferromagnets namely CoFeB and NiFe. For two-sublattice ferromagnet, the inter-sublattice exchange energies become important. Here, we describe the inertial effects in a two-sublattice ferromagnet coupled by the Heisenberg exchange interaction. We follow the similar procedure of Ref.~\cite{Mondal2020nutation} and derive the magnetic susceptibility. We {\it not only} consider the intra-sublattice inertial dynamics, {\it but also} the inter-sublattice dynamics. Our results suggest that there are two precession resonance peaks: one at GHz regime and another at THz regime. Similarly, two nutation peaks can also be observed, both are at the THz regime. By calculating the precession and nutation resonance frequencies, we observe that the resonance frequencies decrease with increasing intra-sublattice relaxation time, however, the scenario is different for inter-sublattice inertial dynamics.

\section{Theory of intra- and inter-sublattice inertial dynamics in two-sublattice ferromagnets}
The inertial dynamics for antiferromagnets have been introduced in Ref.~\cite{Mondal2020nutation}.
For two-sublattice magnetic systems having magnetization  ${\bm M}_{A}$ and  ${\bm M}_{B}$, for $A$ and $B$ representing the two-sublattice, the ILLG equations of motion can be recast as
\begin{align}
\label{sublatticeA}
    \frac{\partial{\bm{M}}_{A}}{\partial t} & = -\gamma_A \left(\bm{M}_A \times \bm{H}_A\right) + \frac{\alpha_{AA}}{M_{A0}}\left(\bm{M}_A\times \frac{\partial{\bm{M}}_{A}}{\partial t}\right) + \frac{\alpha_{AB}}{M_{B0}}\left(\bm{M}_A\times \frac{\partial{\bm{M}}_{B}}{\partial t}\right) \nonumber\\
    & + \frac{\eta_{AA}}{M_{A0}}\left(\bm{M}_A\times \frac{\partial^2{\bm{M}}_{A}}{\partial t^2}\right)  + \frac{\eta_{AB}}{M_{B0}}\left(\bm{M}_A\times \frac{\partial^2{\bm{M}}_{B}}{\partial t^2}\right)\\
    \frac{\partial{\bm{M}}_{B}}{\partial t} & = -\gamma_B \left(\bm{M}_B \times \bm{H}_B\right) + \frac{\alpha_{BB}}{M_{B0}}\left(\bm{M}_B\times \frac{\partial{\bm{M}}_{B}}{\partial t}\right)  + \frac{\alpha_{BA}}{M_{A0}}\left(\bm{M}_B\times \frac{\partial{\bm{M}}_{A}}{\partial t}\right) \nonumber\\
    &+ \frac{\eta_{BB}}{M_{B0}}\left(\bm{M}_B\times \frac{\partial^2{\bm{M}}_{B}}{\partial t^2}\right) + \frac{\eta_{BA}}{M_{A0}}\left(\bm{M}_B\times \frac{\partial^2{\bm{M}}_{A}}{\partial t^2}\right)
    \label{sublatticeB}
\end{align}
In each ILLG dynamics, the first term represents the spin precession around an effective field ${\bm H}_{A/B}$. The intra- and inter-sublattice Gilbert damping dynamics have been denoted by the second and third terms, respectively. Similarly, the last two terms define inertial dynamics. While the intra-sublattice Gilbert and inertial dynamics have been weighed by $\alpha_{AA/BB}$ and $\eta_{AA/BB}$, the same for inter-sublattice dynamics are denoted by  $\alpha_{AB/BA}$ and $\eta_{AB/BA}$. From a simple dimension analysis, it is clear to show that the Gilbert damping parameters $\alpha$ are dimensionless, in contrast, the inertial relaxation times $\eta$ have a dimension of time~\cite{Ciornei2011,Mondal2017Nutation}. It is worth mentioning that the Gilbert damping $\alpha$ has been calculated for several materials within {\it ab initio} frameworks~\cite{kambersky70,kambersky76,kunes02,KunesErratum,tserkovnyak02,steiauf05,Fahnle2006,kambersky07,gilmore07,Brataas2008,Gilmore2008JAP,EbertPRL2011,Thonig_2014,Edwards_2016}, while there are also proposals to calculate the inertial relaxation time within extended breathing Fermi surface model~\cite{Fahnle2011,fahnle2013erratum,Thonig2017}.
These ILLG equations have been contemplated to forecast the signatures of inertial dynamics in collinear antiferromagnets and ferrimagnets~\cite{Mondal2020nutation}.

We consider that the two-sublattice ferromagnet is aligned collinear at the ground state such that ${\bm M}_{A} = M_{A0} \hat{\bm{z}}$ and ${\bm M}_{B} = M_{B0} \hat{\bm{z}}$. The ferromagnetic system is under the application of an external Zeeman field ${\bm H}_0 = H_0\hat{\bm{z}}$. Then, the free energy of the considered two-sublattice system can be considered as the sum of Zeeman, anisotropy, and exchange energies as
\begin{align}
    \mathcal{F}\left(\bm{M}_{A}, \bm{M}_{B}\right)  & =  - H_0\left( M_{Az}+ M_{Bz}\right)  - \frac{K_A}{M^{2}_{A0}} M^2_{Az}  - \frac{K_B}{M^{2}_{B0}} M^2_{Bz} - \frac{J}{M_{A0}M_{B0}}  \bm{M}_A\cdot \bm{M}_B\,,
     \label{Free-energy}
\end{align}
where $K_A$ and $K_B$ are anisotropy energies and $J$ is the isotropic Heisenberg exchange with $J>0$ for ferromagnetic coupling.
To calculate the linear response properties of the system, we consider that the small deviations of magnetization $\bm{m}_{A}(t)$ and $\bm{m}_{B}(t)$ with respect to the ground state are induced by the transverse external field $\bm{h}_{A}(t)$ and $\bm{h}_{B}(t)$. We calculate the effective field in the ILLG equation as the derivative of free energy in Eq.\ (\ref{Free-energy}) to the corresponding magnetization
\begin{align}
\label{HeffA}
     \bm{H}_A  = -\frac{\partial \mathcal{F}\left(\bm{M}_A, \bm{M}_B\right)}{\partial \bm{M}_A} & = \left( H_0 + \frac{2K_A}{M^{2}_{A0}}M_{Az}\right)\hat{\bm{z}}  + \frac{J}{M_{A0}M_{B0}}\bm{M}_B\nonumber\\
     & = \frac{1}{M_{A0}}\left(H_0M_{A0} + 2K_A + J\right) \hat{\bm{z}} + \frac{J}{M_{A0}M_{B0}}\bm{m}_B\,,\\
     \bm{H}_B  = -\frac{\partial \mathcal{F}\left(\bm{M}_A, \bm{M}_B\right)}{\partial \bm{M}_B}  & = \left( H_0 + \frac{2K_B}{M^{2}_{B0}}M_{Bz}\right)\hat{\bm{z}} + \frac{J}{M_{A0}M_{B0}}\bm{M}_A\nonumber\\
     & = \frac{1}{M_{B0}}\left(H_0M_{B0} + 2K_B + J\right) \hat{\bm{z}} + \frac{J}{M_{A0}M_{B0}}\bm{m}_A\,.
     \label{HeffB}
\end{align}
We then expand the magnetization around the ground state in small deviations, $\bm{M}_{A} = M_{A0}\hat{\bm{z}}+\bm{m}_{A} (t)$ and $\bm{M}_{B} = M_{B0}\hat{\bm{z}}+\bm{m}_{B}(t)$. Essentially, with the effective fields in Eqs.~(\ref{HeffA}) and ~(\ref{HeffB}) along with the magnetization, the linear response for sublattice A provides
\begin{align}
  \frac{\partial{\bm{m}}_{A}}{\partial t}
    & = -\frac{\gamma_A}{M_{A0}} \left(H_0M_{A0} + 2K_A + J\right)\left[ m_{Ay} \hat{\bm{x}} - m_{Ax} \hat{\bm{y}} \right] - \frac{\gamma_A J}{M_{B0}}\left[ m_{Bx} \hat{\bm{y}} - m_{By} \hat{\bm{x}}\right] \nonumber\\
    & - \gamma_A M_{A0} \left[ h_{Ax}\hat{\bm{y}} - h_{Ay} \hat{\bm{x}} \right] + \alpha_{AA}\left[\frac{\partial m_{Ax}}{\partial t}\hat{\bm{y}} - \frac{\partial m_{Ay}}{\partial t}\hat{\bm{x}}  \right] + \frac{\alpha_{AB}M_{A0}}{M_{B0}}\left[\frac{\partial m_{Bx}}{\partial t}\hat{\bm{y}} - \frac{\partial m_{By}}{\partial t}\hat{\bm{x}}  \right]\nonumber\\
    & + \eta_{AA} \left[\frac{\partial^2 m_{Ax}}{\partial t^2}\hat{\bm{y}} - \frac{\partial^2 m_{Ay}}{\partial t^2}\hat{\bm{x}}  \right] + \frac{\eta_{AB}M_{A0}}{M_{B0}}\left[\frac{\partial^2 m_{Bx}}{\partial t^2}\hat{\bm{y}} - \frac{\partial^2 m_{By}}{\partial t^2}\hat{\bm{x}}  \right]\,,
\end{align}
obtaining the dynamics for two components $x$ and $y$ as
\begin{align}
    \gamma_AM_{A0} h_{Ax} & = \frac{\gamma_A}{M_{A0}} \left(H_0M_{A0} + 2K_A + J\right) m_{Ax}   - \frac{\gamma_A J}{M_{B0}} m_{Bx}
      + \alpha_{AA}\frac{\partial m_{Ax}}{\partial t} + \frac{\alpha_{AB}M_{A0}}{M_{B0}}\frac{\partial m_{Bx}}{\partial t} \nonumber\\
    &  - \frac{\partial m_{Ay}}{\partial t}+ \eta_{AA} \frac{\partial^2 m_{Ax}}{\partial t^2}  + \frac{\eta_{AB}M_{A0}}{M_{B0}}\frac{\partial^2 m_{Bx}}{\partial t^2}\,,\\
    \gamma_AM_{A0} h_{Ay} & = \frac{\gamma_A}{M_{A0}} \left(H_0M_{A0} + 2K_A + J\right) m_{Ay}   - \frac{\gamma_A J}{M_{B0}} m_{By}
      + \alpha_{AA}\frac{\partial m_{Ay}}{\partial t} + \frac{\alpha_{AB}M_{A0}}{M_{B0}}\frac{\partial m_{By}}{\partial t} \nonumber\\
    &  + \frac{\partial m_{Ax}}{\partial t}+ \eta_{AA} \frac{\partial^2 m_{Ay}}{\partial t^2}  + \frac{\eta_{AB}M_{A0}}{M_{B0}}\frac{\partial^2 m_{By}}{\partial t^2}\,.
\end{align}
In the circular basis defined by $m_{A\pm} = m_{Ax} \pm {\rm i}m_{Ay}$ and $h_{A\pm} = h_{Ax} \pm {\rm i}h_{Ay}$, the equations can be put together
\begin{align}
    \gamma_AM_{A0} h_{A\pm} & = \frac{\gamma_A}{M_{A0}} \left(H_0M_{A0} + 2K_A + J\right) m_{A\pm}   - \frac{\gamma_A J}{M_{B0}} m_{B\pm}
      + \alpha_{AA}\frac{\partial m_{A\pm}}{\partial t} + \frac{\alpha_{AB}M_{A0}}{M_{B0}}\frac{\partial m_{B\pm}}{\partial t} \nonumber\\
    &  \pm {\rm i}\frac{\partial m_{A\mp}}{\partial t}+ \eta_{AA} \frac{\partial^2 m_{A\pm}}{\partial t^2}  + \frac{\eta_{AB}M_{A0}}{M_{B0}}\frac{\partial^2 m_{B\pm}}{\partial t^2}\,.
    \label{h_A}
\end{align}
Similarly, one can calculate the linear response of the sublattice B in the circular basis defined by $m_{B\pm} = m_{Bx} \pm {\rm i}m_{By}$ and $h_{B\pm} = h_{Bx} \pm {\rm i}h_{By}$ as
\begin{align}
    \gamma_B M_{B0} h_{B\pm} & = \frac{\gamma_B}{M_{B0}} \left(H_0M_{B0} + 2K_B + J\right) m_{B\pm}   - \frac{\gamma_B J}{M_{A0}} m_{A\pm}
      + \alpha_{BB}\frac{\partial m_{B\pm}}{\partial t} + \frac{\alpha_{BA}M_{B0}}{M_{A0}}\frac{\partial m_{A\pm}}{\partial t} \nonumber\\
    &  \pm {\rm i}\frac{\partial m_{B\mp}}{\partial t}+ \eta_{BB} \frac{\partial^2 m_{B\pm}}{\partial t^2}  + \frac{\eta_{BA}M_{B0}}{M_{A0}}\frac{\partial^2 m_{A\pm}}{\partial t^2}\,.
    \label{h_B}
\end{align}
We define the response functions $m_{A\pm},m_{B\pm},h_{A\pm},h_{B\pm}  \propto e^{\pm {\rm i} \omega t}$ and  $\Omega_A = \frac{\gamma_A}{M_{A0}} \left(H_0M_{A0} + 2K_A + J\right)$ and $\Omega_B = \frac{\gamma_B}{M_{B0}} \left(H_0M_{B0} + 2K_B + J\right)$.
To simplify the expressions, we introduce the following: $\Gamma_{AA} = \gamma_A M_{A0}$, $\Gamma_{BB} = \gamma_B M_{B0}$, $\Gamma_{AB} = \gamma_A M_{B0}$ and $\Gamma_{BA} = \gamma_B M_{A0}$ such that $\Gamma_{AA}\Gamma_{BB} =\Gamma_{AB}\Gamma_{BA} $. The linear response Eqs.~(\ref{h_A}) and (\ref{h_B}) can be written in a matrix formalism
\begin{align}
    & \begin{pmatrix}
      h_{A\pm}\\
       h_{B\pm}
    \end{pmatrix} \nonumber\\
    & = \begin{pmatrix}
      \dfrac{1}{ \Gamma_{AA}}\left(\Omega_A \pm {\rm i}\omega \alpha_{AA} - \omega^2 \eta_{AA} - \omega\right) & -\dfrac{1}{\Gamma_{AB}} \left(\dfrac{\gamma_A J}{M_{A0}}
       \mp {\rm i}\omega \alpha_{AB}    + \omega^2 \eta_{AB}\right)\\
     -\dfrac{1}{\Gamma_{BA}} \left(\dfrac{\gamma_B J}{M_{B0}}
       \mp {\rm i}\omega \alpha_{BA}    + \omega^2 \eta_{BA}\right) &  \dfrac{1}{ \Gamma_{BB}}\left(\Omega_B \pm {\rm i}\omega \alpha_{BB} - \omega^2 \eta_{BB} - \omega\right)
    \end{pmatrix}\begin{pmatrix}
      m_{A\pm}\\
      m_{B\pm}
    \end{pmatrix}\,.
\end{align}
For finding the susceptibility, we recall $\bm{m}_{\pm} = \chi_{\pm}\cdot {\bm h}_{\pm}$ such that the susceptibility matrix derives as
\begin{align}
   \chi^{AB}_\pm & = \dfrac{1}{\mathcal{D}_\pm}\begin{pmatrix}
      \dfrac{1}{ \Gamma_{BB}}\left(\Omega_B \pm {\rm i}\omega \alpha_{BB} - \omega^2 \eta_{BB} - \omega\right)  & \dfrac{1}{\Gamma_{BA}} \left(\dfrac{\gamma_B J}{M_{B0}}
       \mp {\rm i}\omega \alpha_{BA}    + \omega^2 \eta_{BA}\right)\\
    \dfrac{1}{\Gamma_{AB}} \left(\dfrac{\gamma_A J}{M_{A0}}
       \mp {\rm i}\omega \alpha_{AB}    + \omega^2 \eta_{AB}\right)&  \dfrac{1}{ \Gamma_{AA}}\left(\Omega_A \pm {\rm i}\omega \alpha_{AA} - \omega^2 \eta_{AA} - \omega\right)
    \end{pmatrix}\,,
\end{align}
where the determinant is expressed as \begin{align}
   \mathcal{D}_\pm & = \frac{1}{\Gamma_{AA}\Gamma_{BB}} \left(\Omega_A \pm {\rm i}\omega \alpha_{AA} - \omega^2 \eta_{AA} - \omega\right)\left(\Omega_B \pm {\rm i}\omega \alpha_{BB} - \omega^2 \eta_{BB} - \omega\right)  \nonumber\\
   &  - \frac{1}{\Gamma_{AB}\Gamma_{BA}} \left(\frac{\gamma_A J}{M_{A0}} \mp {\rm i}\omega \alpha_{AB}    + \omega^2 \eta_{AB}\right)\left(\frac{\gamma_B J}{M_{B0}} \mp {\rm i}\omega \alpha_{BA}    + \omega^2 \eta_{BA}\right)\,.
\end{align}

Note that the intra-sublattice dynamical parameters enter in the diagonal elements of the susceptibility matrix, however, the inter-sublattice dynamics are reflected in the off-diagonal elements. Such a susceptibility matrix has been obtained with intra- and inter-sublattice Gilbert damping dynamics for antiferromagnets~\cite{Kamra2018}.
To find the resonance frequencies, one has to solve the equation setting $\mathcal{D}_{\pm} = 0$.
Therefore, a fourth-order equation in frequency is obtained
\begin{align}
    \mathbb{A}_\pm \omega^4 + \mathbb{B}_\pm \omega^3 + \mathbb{C}_\pm \omega^2 + \mathbb{D}_\pm \omega +
    \mathbb{E}_\pm & = 0\,,
    \label{4th-order-Eq}
\end{align}
with the following coefficients
\begin{align}
    \mathbb{A}_\pm & = \eta_{AA}\eta_{BB} - \eta_
    {AB}\eta_{BA}\,,\\
    \mathbb{B}_\pm & =(\eta_{AA} + \eta_{BB}) \mp {\rm i} \left(\alpha_{AA}\eta_{BB} + \alpha_{BB}\eta_{AA}\right) \pm {\rm i}\left(\alpha_{AB}\eta_{BA}+ \alpha_{BA}\eta_{AB}\right)\,,\\
    \mathbb{C}_\pm & =1 \mp {\rm i} \left(\alpha_{AA} + \alpha_{BB}\right)- \left(\Omega_A\eta_{BB} + \Omega_B\eta_{AA}\right) - \alpha_{AA}\alpha_{BB}\nonumber\\
    & -\left(\frac{\gamma_A}{M_{A0}}\eta_{BA} + \frac{\gamma_B}{M_{B0}}\eta_{AB}\right)J-\alpha_{AB}\alpha_{BA}\,,\\
    \mathbb{D}_\pm & =-\left(\Omega_A + \Omega_B\right)\pm {\rm i} \left(\Omega_A\alpha_{BB} +\Omega_B\alpha_{AA} \right) \pm{\rm i} \left(\frac{\gamma_A}{M_{A0}}\alpha_{BA} + \frac{\gamma_B}{M_{B0}}\alpha_{AB}\right)J\,,\\
    \mathbb{E}_\pm & =\Omega_A \Omega_B -\frac{\gamma_A\gamma_B}{M_{A0}M_{B0}}J^2\,.
\end{align}
The analytical solution of the above-mentioned equation is very cumbersome. Therefore, we adopt the numerical techniques for solving Eq.~(\ref{4th-order-Eq}).
The solution of the above equation results in four frequencies, two of them correspond to the precession resonance ($\omega_{\rm p}$) of each sublattice and the other two belong to the nutation resonance ($\omega_{\rm n}$). The real and imaginary parts of the resonance frequency are denoted by ${\tt Re}$ and ${\tt Im}$, respectively. For example, the precession resonance frequencies are $\omega_{\rm p} = {\tt Re}\left(\omega_{\rm p}\right) + {\rm i}{\tt Im}\left(\omega_{\rm p}\right)$, while the nutation resonance frequencies are $\omega_{\rm n} = {\tt Re}\left(\omega_{\rm n}\right) + {\rm i}{\tt Im}\left(\omega_{\rm n}\right)$.
Comparing Eq.~(\ref{4th-order-Eq}), a similar equation has been obtained for antiferromagnets and ferrimagnets~\cite{Mondal2020nutation}, however, without the inter-sublattice inertial dynamics. We mention that the inter-sublattice Gilbert damping dynamics have extensively been discussed~\cite{Kamra2018,Yuan_2019}. Therefore, we will not consider in the following discussions. In particular, we allow $\alpha_{AB} = \alpha_{BA} = 0$, and calculate the inertial effects on precession and nutation resonances.

\section{Numerical results}
To calculate the resonace frequencies, we numerically solve the Eq.~(\ref{4th-order-Eq}) for two-sublattice ferromagnets having same magnetic moments in each sublattice i.e., $M_{A0} = M_{B0}$.
We use the following parameters: $\gamma_A = \gamma_B = 1.76\times 10^{11}$ T$^{-1}$-s$^{-1}$, $J = 10^{-21}$ J, $K_A = K_B = 10^{-23}$ J, $\alpha_{AA} = \alpha_{BB} =  \alpha =  0.05$, $\alpha_{AB} = \alpha_{BA} = 0$. The considered exchange and anisotropy energies have similar order of magnitude as typical ferromagnets e.g., Fe \cite{Pajda2001}. The chosen Gilbert damping $\alpha = 0.05$ is within the {\it ab initio} reported values~\cite{EbertPRL2011}. For inertial relaxation times, even though, the {\it ab initio} calculation suggests about fs timescales for transition metals~\cite{Thonig2017}, the recent experiment predicts it to be a higher value up to several hundreds of fs~\cite{neeraj2019experimental}. Therefore, in what follows, we have considered the inertial relaxation times ranging from fs to ps.

\subsection{Intra-sublattice inertial dynamics}
\begin{figure}[hbt!]
    \centering
    \includegraphics[scale = 0.22]{Fig2.pdf}
    \caption{The calculated resonance frequencies  as a function of intra-sublattice inertial relaxation time for two-sublattice ferromagnets using $M_{A0} = M_{B0} = 2\mu_B$. (a) The precession and nutation resonance frequencies and (b) the effective Gilbert damping have been plotted.}
    \label{intra-sublattice-Fig1}
\end{figure}
\begin{figure}[hbt!]
    \centering
    \includegraphics[scale = 0.9]{Fig3.pdf}

    \caption{The calculated precession resonance frequencies as a function of $M_{A0}/M_{B0}$ for two-sublattice ferromagnets, at several intra-sublattice inertial relaxation times. (a) The real part of the lower precession resonance frequencies,  (b) the effective damping of lower resonance mode, (c) the real part of the upper precession resonance frequencies,  (d) the effective damping of upper resonance mode, has been plotted.}
    \label{intra-sublattice-Fig2}
\end{figure}
 To focus on the intra-sublattice inertial dynamics, we set the inter-sublattice relaxation time to zero i.e., $\eta_{AB} = \eta_{BA}$ = 0, keeping the same inertial relaxation time in two-sublattice $\eta_{AA} = \eta_{BB} = \eta$. With this set of specifications, the calculated frequencies have been shown in Fig.~\ref{intra-sublattice-Fig1}. One can see that there exist two precession resonance frequencies (positive) and the corresponding two nutation resonance frequencies (negative). We denote these two positive precession frequencies as $\omega_{\rm p+}^{\rm u}$ and $\omega_{\rm p+}^{\rm l}$, while the two negative nutation frequencies are $\omega_{\rm n-}^{\rm u}$ and $\omega_{\rm n-}^{\rm l}$. The superscripts ``u'' and ``l'' denote the upper and lower frequencies, respectively. These results are in contrast with the observation in antiferromagnets or ferrimagnets, where one positive and one negative precession (and nutation) frequencies are expected \cite{Mondal2020nutation}. { Nevertheless, the quantitative comparison of the calculated frequencies agrees with those of the ferrimagnets, where the upper (THz), and lower (GHz) frequency precession resonances are called an exchange and ferromagnetic modes, respectively~\cite{Mondal2020nutation,Schlickeiser2012}.  Similar to antiferromagnets and ferrimagnets~\cite{Mondal2020nutation}, the resonance frequencies decrease with the intra-sublattice inertial relaxation time in the case of two-sublattice ferromagnets. Especially, the lower nutation resonance frequency scales with $1/\eta$, while the upper one shows deviation from $1/\eta$ at higher relaxation times. This deviation from $1/\eta$ has been noticed in two nutation modes for antiferromagnets and ferrimagnets~\cite{Mondal2020nutation}. An interesting feature is that the precession and nutation frequencies cross each other at certain inertial relaxation times in ferromagnets. Such crossing was not observed in antiferromagnets and ferrimagnets~\cite{Mondal2020nutation}. The crossing happens especially with the upper precession mode with lower nutation mode as seen in Fig.~\ref{intra-sublattice-Fig1}(a). However, we note that crossing of these two modes have positive and negative frequencies, meaning that the upper precession mode ($\omega_{\rm p+}^{\rm u}$) has a positive rotational sense, however, the lower nutation mode ($\omega_{\rm n-}^{\rm l}$) has the opposite rotational sense in circular basis.

The inertial dynamics affect the effective Gilbert damping in a system. This has been demonstrated in Fig.~\ref{intra-sublattice-Fig1}(b) for two-sublattice ferromagnet by the ratio of imaginary and real parts of the calculated frequencies. We have used the same Gilbert damping for both the sublattices $\alpha \sim 0.05$ and therefore, the effective damping remains the same at smaller inertial relaxation times. However, the effective damping decreases with increased relaxation times, a fact that is consistent with the results of antiferromagnets \cite{Mondal2020nutation}. It is observed that the decrease in effective damping is exactly the same for precession and corresponding nutation modes. Moreover, the upper precession mode is influenced strongly, which has already been observed for ferrimagnets \cite{Mondal2020nutation}.

Next, we calculate the influence of different sublattice magnetic moment ($M_{A0} \neq M_{B0}$) on inertial dynamics. In particular, we compute the precession resonance frequencies as a function of the ratio of magnetic moments ($M_{A0}/M_{B0}$), at several inertial relaxation times in Fig.~\ref{intra-sublattice-Fig2}. We observe that the resonance frequencies decrease with increasing difference in the magnetic moments. Such reduction is less visible in case of lower precession frequencies e.g., Fig.~\ref{intra-sublattice-Fig2}(a), however, more prominent in upper precession frequencies in Fig.~\ref{intra-sublattice-Fig2}(c). However, the difference of frequencies calculated at several relaxation times are similar for $M_{A0} = M_{B0}$ and $M_{A0} \neq M_{B0}$.
The latter suggests that the inertial dynamics do not get quantitatively influenced by the same or different sublattice magnetic moments. A similar conclusion can also be made from the computation of effective damping in Figs.~\ref{intra-sublattice-Fig2}(b) and \ref{intra-sublattice-Fig2}(d). The effective damping for the upper and lower precession modes remains almost constant (with a very small positive slope) for a higher ratio of $M_{A0} / M_{B0}$.

\subsection{Inter-sublattice inertial dynamics}
\begin{figure}[hbt!]
    \centering
    \includegraphics[scale = 0.22]{Fig4.pdf}
    \caption{The calculated resonance frequencies  as a function of inter-sublattice inertial relaxation time for two-sublattice ferromagnets using $M_{A0} = M_{B0} = 2\mu_B$. The intra-sublattice inertial relaxation time was kept constant $\eta = 100$ fs.  (a) The precession and nutation resonance frequencies and (b) the effective Gilbert damping have been plotted.}
    \label{inter-sublattice-Fig3}
\end{figure}
To investigate the inter-sublattice inertial dynamics, we set the intra-sublattice relaxation time as $\eta_{AA} = \eta_{BB} = \eta = 100$ fs. Such a relaxation time is lower than the experimental findings in two-sublattice ferromagnets~\cite{neeraj2019experimental}. In fact, the direct comparison of Eq.~(\ref{sublatticeA}) with the Eq.~(2) of Ref.~\cite{neeraj2019experimental} provides $\eta \sim \alpha \tau$. With the experimental findings for CoFeB, $\alpha = 0.0044$ and $\tau = 72$ ps (see Table 1 in Ref.~\cite{neeraj2019experimental}), we calculate $\eta = 316$ fs. We compute the effect of inter-sublattice inertial dynamics as a function of $\eta_{AB} = \eta_{BA} = \eta^\prime$ in Fig.~\ref{inter-sublattice-Fig3} considering $\eta^\prime < \eta$.
As we mentioned earlier, the overlapping of precession ($\omega_{\rm p+}^{\rm u}$) and nutation ($\omega_{\rm n-}^{\rm l}$) frequencies at the intra-sublattice relaxation time $\eta = 100$ fs can be seen.
We observe that the upper precession resonance frequency ($\omega_{\rm p+}^{\rm u}$) increases, while the lower one ($\omega_{\rm p+}^{\rm l}$) decreases very small with inter-sublattice relaxation times. A similar conclusion can be made for nutation frequencies. This is in contrast to the observation of intra-sublattice inertial dynamics as discussed above. A divergence in the upper nutation frequency can be noticed at the limit $\eta^\prime \rightarrow \eta$. Such divergence can be explained through the coefficient $\mathbb{A}$ in Eq.~(\ref{4th-order-Eq}). At the limit $\eta^\prime \rightarrow \eta$, the coefficient of fourth power in frequency becomes $\mathbb{A} = \eta_{AA}\eta_{BB} - \eta_{AB}\eta_{BA} = \eta^2 - \eta^{\prime 2} \rightarrow 0$, which brings the fourth-order equation into an effective third-order equation in frequency.

A similar observation can also be concluded from the calculation of effective damping in Fig.~\ref{inter-sublattice-Fig3}(b). Similar to the intra-sublattice inertial dynamics, the effective damping of the precession and corresponding nutation mode behaves exactly the same for the inter-sublattice inertial dynamics. We observe that the damping of upper precession and nutation modes increases with inter-sublattice inertial relaxation time, however, it is the opposite for lower precession and nutation modes. Therefore, we conclude that the effect of intra- and inter-sublattice inertial dynamics are contrasting.

\section{Conclusions}
To conclude, we have incorporated the intra- and inter-sublattice inertial dynamics within the LLG equation of motion and calculated the FMR resonance for two-sublattice ferromagnets. To this end, we first derive the magnetic susceptibility that is a tensor. To calculate the resonance frequencies, we find the poles of the susceptibility. Without the inertial dynamics, there exist two precession modes in a typical two-sublattice ferromagnet. The introduction of inertial dynamics shows two nutation resonance frequencies corresponding to the precession modes. We note that these precession and nutation resonances can be excited by right and left circularly polarised pulses, respectively, and vice-versa within a circular basis. The precession and nutation frequencies decrease with the intra-sublattice relaxation time as also has been seen in the case of antiferromagnets in previous work~\cite{Mondal2020nutation}. However, at certain relaxation times, the precession and nutation frequencies overlap with each other. Note that these overlapping precession and nutation frequencies have opposite rotational sense in circular basis, thus, they can be neatly realised in the experiments. The inter-sublattice inertial dynamics increase the resonance frequencies and effective damping for upper precession mode, however, have opposite effect on lower precession mode in two-sublattice ferromagnets.

\section{Acknowledgments}
The author acknowledges Levente R\'ozsa and Ulrich Nowak for valuable discussions and the Swedish Research Council (VR 2019-06313) for research funding.

\bibliographystyle{iopart-num}

\providecommand{\newblock}{}