\subsection{SINR of each link}
Knowing the distance $r_k(t)$ between the two communicating nodes, the SINR of each link in LoS/NLoS can be defined as
\begin{equation}
    \chi_{k,n|r_{k}(t)}(t) = \frac{P_k G_{k,M} L^{-1}_{k,n}(r_k(t))|h_{k,n|r_k(t)}|^2}{I_{k}(t) + \sigma^2}, \tag{9}
\end{equation}
with $P_{SD} = P_{SR} = P_T$, $P_{RD} = P_R$, $h_{k,n|r_k(t)}$ being the small-scale fading of the considered link, and
\begin{align}
   & L_{SD,n}(r) = L_{Gn}(1+r), \nonumber\\
   & L_{SR,n}(r) = L_{RD,n}(r) = L_{An}\left(\sqrt{r^2+H_R^2}\right). \nonumber
\end{align}
The aggregate interferences of the three links are given by
\begin{equation}
I_k(t) = \sum_{\bm{\phi}_j \in \bm{\Phi}_{k,I}(t)} P_k G_{k,m} L_{k,n_j}^{-1}(\rVert \bm{\phi}_j \rVert)|h_{k,n_j|\rVert\bm{\phi}_j\rVert}|^2, \tag{10}
\end{equation}
with $n_j = \{L,N\}$ depending on whether the link between $\bm{\phi}_j$ and the node of interest is in LoS or in NLoS, and $h_{k,n_j,j|\rVert\bm{\phi}_j\rVert}$ being the small-scale fading of this link.\smallskip

The interferers' PPs are defined as
\begin{equation}
\bm{\Phi}_{k,I}(t) = \bm{\Phi}_k \backslash \{\tilde{\bm{\phi}}_{k}(t)\},\tag{11}
\end{equation}
with $\bm{\Phi}_{SD} = \bm{\Phi}_{SR} = \bm{\Phi}_T$, and $\bm{\Phi}_{RD} = \bm{\Phi}_R$. Even if the selected TBS can vary at different time instants, with instantaneous handovers, \cite{Banagar2019} has shown that the distribution of $r_{SR}(t)$ is time-invariant.\smallskip

Knowing $r_{SD}$ and $r_{SR}(t)$, the densities of the interferers' PPs for the direct link and the first hop of the relay link are given by
\begin{alignat}{10}
&\lambda_{SD,I|r_{SD}}(r,t) &&=\lambda_T \: u(r-r_{SD}), \tag{12}\\
&\lambda_{SR,I|r_{SR}(t)}(r,t) &&= \lambda_T \: u(r-r_{SR}(t)). \tag{13}
\end{alignat}

\begin{figure*}[b]
\rule{\linewidth}{0.5pt}
\begin{multline}
  \mathcal{A}_{SRD} = \exp{-\pi\lambda_T H_R^2}\left(1-\exp{-\pi\lambda_R\: v^2t^2}\right) \\+ 2\pi\lambda_R\:\exp{-\pi\lambda_T(H_R^2+v^2t^2)}\int_{vt}^\infty r\:\exp{-\pi(\lambda_R+\lambda_T)r^2 + 2\pi\lambda_T vt\:r}\mathrm{d}r\label{eqn:probsel} \tag{18}
\end{multline}
\end{figure*}
\vspace*{-0.0cm}

\subsection{PDF of the initial distances between the selected nodes}
Since the TBSs' and UAVs' HPPPs are independent, $r_{SD}$ and $r_{RD}(t)$ are independent. It is also the case for $r_{SR}(t)$ and $r_{RD}(t)$. However, $r_{SR}(t)$ is dependent of $r_{SD}$ since there is a non zero probability for the same TBS to be selected by the typical UE and the selected RN, as illustrated by Lemma 1.
\medskip

\textit{\textbf{Lemma 1:} The probability for the same TBS to be selected by the typical UE and the selected RN, conditioned on $r_{SD}$ and $r_{RD,0}$ is given by
\begin{align}
    &\Pr{\bm{\tilde{\phi}}_{SD},\bm{\tilde{\phi}}_{SR}(t)|r_{SD},r_{RD,0}}{\bm{\tilde{\phi}}_{SD} = \bm{\tilde{\phi}}_{SR}(t)}\label{eqn:probOfdependency} \tag{14} \\
    &= \exp{-\pi\lambda_T\left(r_{SD}^2+r_{RD}^2(t)\right)}I_0\left(2\pi\lambda_T\:r_{SD}r_{RD}(t)\right). \nonumber
\end{align}}
\textit{Proof:} see Appendix A. \hfill{\small$\blacksquare$}\medskip

Therefore, in order to simplify the mathematical expressions, Approximation 1 is introduced.\medskip

\textit{\textbf{Approximation 1:} We assume that $r_{SR}(t)$ is independent with $r_{SD}$.}
\medskip

Since \eqref{eqn:probOfdependency} decreases with $\lambda_T$, Approximation 1 becomes more accurate when the TBS' density increases because the dependence between $\tilde{\bm{\phi}}_{SD}$ and $\tilde{\bm{\phi}}_{SR}(t)$ decreases.\medskip

\textit{\textbf{Lemma 2:} With the nearest neighbour association rule, the PDFs of the initial distances between the couple of nodes of the first hop, second hop and direct link are given by
\begin{equation}
   f_{r_{k,0}}(r_{k,0}) = 2\pi\lambda_k\:r_{k,0}\: \exp{-\pi\lambda_k r_{k,0}^2},\label{eqn:PDFrSD} \tag{15}
\end{equation}
with $\lambda_{SD} = \lambda_{SR} = \lambda_T,\ \lambda_{RD} = \lambda_R$.}\smallskip

\textit{Proof:} see Appendix B. \hfill{\small$\blacksquare$}\medskip

With the nearest neighbour association rule, we define the regions of association $\Omega_{SD}(t)$ and $\Omega_{SRD}(t)$ as
\begin{alignat}{10}
    &\Omega_{SD}(t) &&:= \left\{(r_{SD},r_{RD,0})\:|\: r_{SD}^2 \leq r_{RD}^2(t) + H_R^2\right\}, \tag{16}\\
    &\Omega_{SRD}(t) &&:= \mathbb{R}_+^2\:\backslash\:\Omega_{SD}(t). \tag{17}
\end{alignat}
Therefore, at time $t\geq 0$, for a given couple of initial distances $\bm{x}=(r_{SD},r_{RD,0})$, the typical UE selects the direct or relay link respectively if $\bm{x}\in\Omega_{SD}(t)$ or $\bm{x}\in\Omega_{SRD}(t)$. This is illustrated in Figure \ref{fig:associationRegionsNNAMob}. The integration of the PDFs on these regions gives the probability of association with the direct and relay link. The result is presented in Lemma 3.\medskip

\begin{figure*}[b]
\rule{\linewidth}{0.5pt}
\begin{alignat}{2}
&\mathcal{I}_{SD}^{(a)}(\beta) &&= \int_0^{H_R}\mathcal{P}_{SD|r_{SD}=r}(\beta)f_{r_{SD}}(r) \:\mathrm{d}r\label{eqn:I1}\tag{20}\\
&\mathcal{I}_{SD}^{(b)}(\beta,t) &&= \exp{-\pi\lambda_R \left(v^2t^2-H_R^2\right)} \int_{H_R}^\infty \mathcal{P}_{SD|r_{SD}=r}(\beta) f_{r_{SD}}(r)\: \exp{-\pi\lambda_R \left(r^2 + 2vt\sqrt{r^2-H_R^2}\right)}\mathrm{d}r\label{eqn:I2}\tag{21}\\
&\mathcal{I}_{SRD}^{(a)}(\beta,t) &&= \exp{-\pi\lambda_T H_R^2}\int_0^{\infty} \mathcal{P}_{SR|r_{SR}(t)=r}(\beta) f_{r_{SR}(t)}(r) \: \mathrm{d}r\int_{0}^{vt} \mathcal{P}_{RD|r_{RD,0}=r}(\beta,t) f_{r_{RD,0}}(r)\:\mathrm{d}r\label{eqn:I3}\tag{22}\\
\nonumber&\mathcal{I}_{SRD}^{(b)}(\beta,t) &&= \exp{-\pi\lambda_T \left(H_R^2+v^2t^2\right)}\int_0^{\infty} \mathcal{P}_{SR|r_{SR}(t)=r}(\beta) f_{r_{SR}(t)}(r) \: \mathrm{d}r\\
&&&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\cdot\int_{vt}^\infty \mathcal{P}_{RD|r_{RD,0}=r}(\beta,t)f_{r_{RD,0}}(r)\:\exp{-\pi\lambda_T (r^2-2vt\:r)}\mathrm{d}r\label{eqn:I4}\tag{23}
\end{alignat}
\rule{\linewidth}{0.5pt}
\begin{equation}
    \mathcal{P}_{k|n,r_{k}(t)}(\beta,t) = \frac{1}{2} +
    \frac{1}{\pi}\int_0^\infty\frac{1}{\tau}
    \mathrm{Im}\left[\mathcal{L}_{|h_{k,n|r_k(t)}|^2}\left(-j\tau P_kG_{k,M}L_{k,n}^{-1}(r_k(t))\right)
    \mathcal{L}_{I_{k}(t)|r_{k}(t)}(j\tau\beta) \: \exp{-j\tau\beta \sigma^2}\right]\mathrm{d}\tau\label{eqn:PSDRSD} \tag{25}
\end{equation}
\rule{\linewidth}{0.5pt}
\begin{equation}
    \mathcal{L}_{I_k(t)|r_k(t)}(s) = \exp{-2\pi \int_{0}^\infty \left(1 - \sum_{n=\{L,N\}} p_{k,n}(r) \mathcal{L}_{|{h}_{k,n|r}|^2}\left(s\:P_k G_{k,m}L_{k,n}^{-1}(r) \right)\right) \lambda_{k,I|r_k(t)}(r,t)\: r\mathrm{d}r} \label{eqn:LIT} \tag{27}
\end{equation}
\end{figure*}

\textit{\textbf{Lemma 3:} With the nearest neighbour association rule, the probability $\mathcal{A}_{SRD}$ that the typical UE selects the relay link at time $t\geq 0$ is given by \eqref{eqn:probsel}.}
\vspace*{0.1cm}

\textit{Proof:} see Appendix C. \hfill{\small$\blacksquare$}\medskip

\begin{figure}[H]
\begin{center}
\begin{adjustbox}{width = 0.8\linewidth}
\begin{tikzpicture}
    \node (RSD) at (1,4) {$r_{SD}$};
    \node (RRD) at (7,1) {$r_{RD,0}$};
    \node[anchor=north] (vt) at (2.5,1) {$vt$};
    \node[anchor=east] (HR)  at (1,2) {$H_R$};
    \draw[->,thick] (1,1) -- (RSD);
    \draw[->,thick] (1,1) -- (RRD);
    \draw[thick,-] (HR) -- (2.5,2);
    \draw[-,black!20!white,dashed] (vt) -- (2.5,2);
    \draw[thick,-,dashed,black!50!white] (2.5,2) -- (2.5,4.2);
    \draw[thick,-,dashed,black!50!white] (2.5,2) -- (6.5,2);
    \node[fill=white] (SD)  at (5.25,2) {$\Omega_{SD}(t)$};
    \node[black!50!white] (ISDA) at (3,1.5) {$\mathcal{I}_{SD}^{(a)}$};
    \node[black!50!white] (ISDB) at (5.5,3) {$\mathcal{I}_{SD}^{(b)}$};
    \node[fill=white] (SRD) at (2.7,3.5) {$\Omega_{SRD}(t)$};
    \node[black!50!white] (ISRDA) at (1.8,2.8) {$\mathcal{I}_{SRD}^{(a)}$};
    \node[black!50!white] (ISRDB) at (3.8,4.2) {$\mathcal{I}_{SRD}^{(b)}$};
    \draw[thick,samples=100,domain=0:2.9,shift={(2.5,2)},
    postaction={decorate,decoration={text along path,text align={right},raise=1ex,text={{\scriptsize $r_{SD}^2=r_{RD}^2(t)+H_R^2$}{}}}}] plot ({\x}, {sqrt(\x^2+1)-1});
\end{tikzpicture}
\end{adjustbox}
\end{center}
\caption{\small Regions of association with the nearest neighbour association rule for the second mobility scheme at time $t\geq 0$. The typical UE selects the direct or the relay link if the couple $(r_{SD},r_{RD,0})$ is respectively in $\Omega_{SD}(t)$ or $\Omega_{SRD}(t)$.}
\label{fig:associationRegionsNNAMob}
\end{figure}

\subsection{Coverage Probability}
Using the definition of the SINRs of each link, Lemmas 2 and Approximation 1, Proposition 1 gives the mathematical expression of the CP of the typical UE.
\medskip

\textit{\textbf{Proposition 1:} With the second mobility scheme, considering the nearest neighbour association rule and selection between direct and relay link, the approximate CP of the typical UE at time $t\geq 0$ for a given SINR threshold $\beta$ is given by
\begin{equation}
\mathcal{P}(\beta,t) = \mathcal{I}_{SD}^{(a)}(\beta) + \mathcal{I}_{SD}^{(b)}(\beta,t) + \mathcal{I}_{SRD}^{(a)}(\beta,t) + \mathcal{I}_{SRD}^{(b)}(\beta,t), \tag{19}
\end{equation}
where $\mathcal{I}_{SD}^{(a)}$, $\mathcal{I}_{SD}^{(b)}$, $\mathcal{I}_{SRD}^{(a)}$ and $\mathcal{I}_{SRD}^{(b)}$ are respectively given by \eqref{eqn:I1}-\eqref{eqn:I4}. The subscripts $SD$ and $SRD$ correspond to whether the considered couples $(r_{SD},r_{RD,0})$ are included in $\Omega_{SD}(t)$ or $\Omega_{SRD}(t)$, as shown in Figure \ref{fig:associationRegionsNNAMob}. In these expressions, $\mathcal{P}_{k|r_{k}(t)}$ are respectively the CP of each link conditioned on the distance between the two nodes of interest at time $t$.}\smallskip

\textit{Proof:} see Appendix D. \hfill{\small$\blacksquare$}\medskip

Moreover, the CPs of every link taken independently are obtained by removing the conditioning of $\mathcal{P}_{k|r_{k}(t)}$. Using the law of total probability, we have
\begin{equation}
    \mathcal{P}_{k|r_{k}(t)}(\beta,t) = \sum_{n=\{L,N\}}p_{k,n}(r_k(t))\:\mathcal{P}_{k|n,r_k(t)}(\beta,t). \tag{24}
\end{equation}
These expressions are functions of the CPs of each hop if these links are in LoS or NLoS, and respectively conditioned on the distances between the nodes of interest. The Gil-Pelaez theorem enables to compute these quantities based on their characteristic functions \cite{GilPelaez}. It enables to reduce the mathematical complexity of the framework (no multiple numerical derivations of the Laplace Transforms (LT) of the aggregate interferences for each link) using numerical integrations. Based on this theorem, Proposition 2 gives the expressions of $\mathcal{P}_{k|n,r_{k}(t)}$.
\medskip

\textit{\textbf{Proposition 2:} The CPs of each link in LoS or NLoS, conditioned on the distance between the two communicating nodes are given by \eqref{eqn:PSDRSD}. In these expressions, $\mathcal{L}_{|h_{k,n|r_k(t)}|^2}$ are the LTs of the small-scale fading channel power gain of each link, while $\mathcal{L}_{I_{k}(t)|r_{k}(t)}$ are the LTs of the aggregate interferences respectively conditioned on $r_{SD}$, $r_{SR}(t)$ and $r_{RD}(t)$.}\smallskip

\textit{Proof:} see Appendix E. \hfill{\small$\blacksquare$}\medskip

\subsection{LTs of the small-scale fading channel power gains}
\textit{\textbf{Lemma 4:} Under Rician fading, the LTs of the small-scale fading channel power gains of each link are given by
\begin{equation}
\mathcal{L}_{|h_{k,n|r}|^2}(s) = \frac{K_{k,n}(r)+1}{K_{k,n}(r)+1+s}\:\exp{-\frac{K_{k,n}(r)\:s}{K_{k,n}(r)+1+s}}.
\label{eqn:LTRicianFading} \tag{26}
\end{equation}
} \smallskip
\textit{Proof:} see Appendix F. \hfill{\small$\blacksquare$}\medskip

\subsection{LTs of the aggregate interferences}
\textit{\textbf{Proposition 3:} The LTs of the aggregate interferences of each link at time $t \geq 0$, respectively conditioned on $r_{SD}$, $r_{SR}(t)$ and $r_{RD}(t)$ are given by \eqref{eqn:LIT}.}\smallskip

\textit{Proof:} see Appendix G. \hfill{\small$\blacksquare$}\medskip
\subsection*{Contributions}
The contributions of this paper are summarised here below:
\begin{enumerate}
    \item We consider a relay network composed of TBSs and UAV RNs in an infinite area, with positions modelled using two independent HPPPs. The UAV RNs are either hovering or are allowed to move following two mobility schemes presented in \cite{Banagar2019}. The association policy is the following: the typical UE is either served via a direct link (from a TBS), or via a relay link (TBS - UAV RN - UE). The selection of the link is performed using a nearest neighbour association rule.
    \item We derive an approximate expression for the Coverage Probability (CP) of the typical UE, and the CP of each link taken independently. We resort to the Gil-Pelaez theorem \cite{GilPelaez} for an increased tractability.
    \item Different path loss and small-scale fading models are used depending on the considered link, in order to make the scenario more realistic. Every link is either in LoS or NLoS, based on a probability model developed for Air-to-Ground (A2G) and Ground-to-Ground (G2G) cases. The parameters of the fading and path loss of each link depend on its LoS/NLoS and A2G/G2G nature. Regarding the small-scale fading, Rician fading is considered, and the introduction of distance-dependent K factors is proposed.
\end{enumerate}
\smallskip

This paper is organised as follows: the system model is detailed in Section \ref{section:system_model}. Then, the coverage probability is derived in Section \ref{section:coverage_probability}. Finally, numerical results are presented in Section \ref{section:numerical_results}.

\subsection{Network topology}
We consider a typical UE located at $(0,0)$ in $\mathbb{R}^2$. The locations of the UAV RNs and TBSs are modelled as two independent HPPPs $\bm{\Phi}_R$ and $\bm{\Phi}_T$ of densities $\lambda_{R}$ and $\lambda_{T}$. The UAVs are all flying at a fixed altitude $H_R$, while the TBSs and the typical UE are at ground level. The transmit powers of the UAVs and TBSs are respectively given by $P_R$ and $P_T$. The network is assumed to be highly congested, such that all RNs and TBSs interfere with each other while transmitting.

\subsection{Mobility Schemes}
The UAVs are allowed to move following two mobility schemes in $\mathbb{R}^2$ presented in \cite{Banagar2019}, both at a fixed altitude.\smallskip

In the first mobility scheme, all the UAVs move independently at a constant speed $v$ in random directions in $\mathbb{R}^2$. It has been shown in \cite{Banagar2019} that the RNs' HPPP remains an HPPP with the same density at every time $t \geq 0$ with this scheme. Therefore, assuming an instantaneous handover between different RNs, the performance at the typical UE is not affected by the movements of the UAV RNs.\smallskip

In the second mobility scheme, all the UAVs move independently at a constant speed $v$ in random directions in $\mathbb{R}^2$, except the serving UAV, which moves towards the typical UE at position $(0,0)$ and then hovers above it. Knowing the initial distance between the selected RN and the typical UE, the density of the RNs interferers' Point Process (PP) is given by \eqref{eqn:densityMob} in that scenario \cite{Banagar2019}. With this mobility scheme, the distance between the typical UE and the selected RN at time $t\geq 0$ is given by
\begin{equation}
    r_{RD}(t) = (r_{RD,0} - vt)\:u(r_{RD,0} - vt), \tag{2}
\end{equation}
with $r_{RD,0}=r_{RD}(0)$, and $u$ being the Heaviside step function.

\subsection{Association Rule}
For sake of mathematical tractability, the nearest neighbour criterion is used as association rule to select the serving nodes:
\begin{equation}
    \tilde{\bm{\phi}} = \argmin_{\bm{\phi}_i \in \bm{\Phi}} \left\{\left\rVert \bm{\phi}_i - \bm{x} \right\rVert \right\}. \tag{3}
\end{equation}
This rule is also applied to make the selection between the best direct link (TBS to UE) and the best relay link (TBS to UE through a RN). Since the selected links can vary at different time instants with mobile UAVs, we assume instantaneous handovers between different UAVs or between direct and relay links.\smallskip

In summary, the typical UE selects the direct link if the distance with its closest TBS $\tilde{\bm{\phi}}_{SD}$ is smaller than the distance with its closest RN $\tilde{\bm{\phi}}_{RD}(t)$. Otherwise the relay link is selected. In that case, the selected RN is served by its closest TBS $\tilde{\bm{\phi}}_{SR}(t)$. This scenario is illustrated in Figure \ref{fig:scenario}.

\begin{figure}[H]
    \centering
    \includegraphics[width=0.7\linewidth]{Images/scenario.PNG}
    \caption{\small Illustration of the scenario: the typical UE is at $(0,0)$ and selects either the direct or the relay link depending on the distances between the nodes. The RNs and TBSs are respectively represented by blue and red dots. The selected nodes are circled in green.}
    \label{fig:scenario}
\end{figure}
\vspace*{-0.3cm}

\subsection{Relay Transmission Protocol}
In this paper, the UAVs adopt the Decode-and-Forward (DF) protocol, executed over two time instants. At the first time instant, all the TBSs are transmitting to UEs or RNs. Thus the TBSs interfere with each other at that moment. At the second time instant, all the RNs re-transmit the signal to the UEs. Therefore, the RNs are interfering with each other at that instant. The Signal-to-Interference-plus-Noise-Ratio (SINR) of the relay link is given by \cite{Laneman2004}
\begin{equation}
    \chi_{SRD}(t) = \min\left(\chi_{SR}(t),\chi_{RD}(t)\right),
    \label{eq:SINRDF} \tag{4}
\end{equation}
where $\chi_{SR}(t)$ is the SINR between the TBS and the RN, and $\chi_{RD}(t)$ is the SINR between the RN and the typical UE.

\subsection{Beamforming model}
The TBSs and UAVs are equipped with antenna arrays to generate directional beams and perform transmit and receive beamforming. This leads to an SINR increase for the same transmitted power, enabling to compensate for a high path loss. As proposed in \cite{Belbase2018,Biswas2016,Ma2019}, a simple sectored antenna model is introduced:
\begin{equation}
    G(\psi,\theta) = \left\{
    \begin{matrix}
    G_{M} &\mathrm{if} |\psi| \leq \frac{\psi_{bw}}{2}\ \mathrm{and}\ |\theta| \leq \frac{\theta_{bw}}{2}\\
    G_{m} &\mathrm{otherwise}
    \end{matrix}
    \right., \tag{5}
\end{equation}
where $G_{M}$ is the maximum beamforming gain obtained when the elevation angle $\theta$ and the azimuth angle $\psi$ are respectively within the half elevation and azimuth power beamwidth $\theta_{bw}$ and $\psi_{bw}$, and $G_{m}$ is the minimum beamforming gain obtained in the other cases. In practice, it will be assumed that the beam of the target receiver is perfectly aligned with the beam of the transmitter and misaligned with the beams of the interferers, leading thereby to the array gains listed in Table \ref{tab:beamforming_gains}. The maximum and minimum beamforming gains of the TBSs and RNs are respectively denoted by $G_{TM}$, $G_{Tm}$, $G_{RM}$ and $G_{Rm}$.
\begin{table}[H]
\begin{center}
\begin{tabular}{|l||l|l|}
    \hline
         & \multicolumn{1}{c|}{Target} & \multicolumn{1}{c|}{Interference} \\\hline\hline
            {TBS-UE} & $G_{SD,M}=G_{TM}$ & $G_{SD,m}=G_{Tm}$ \\\hline
        {TBS-RN} & $G_{SR,M}=G_{TM}G_{RM}$ & $G_{SR,m}=G_{Tm}G_{Rm}$ \\\hline
        {RN-UE} & $G_{RD,M} = G_{RM}$ & $G_{RD,m} = G_{Rm}$ \\\hline
    \end{tabular}
    \caption{\small Beamforming gains for the different links.}
    \label{tab:beamforming_gains}
\end{center}
\end{table}
\vspace*{-0.6cm}
\subsection{Propagation model}
The channel model takes into account path loss and small-scale fading. The parameters of the different models are different for A2G and G2G links, in LoS or NLoS. We also introduce a model for the LoS probability.\smallskip

\paragraph{Path loss}
We consider two linear path loss models for A2G and G2G links \cite{Khawaja2019}:
\begin{equation}
    L_{ln}(r)\ [dB] = \beta_{ln} + 10\:\alpha_{ln} \: \log_{10}\left(r\right), \tag{6}
\end{equation}
where $r$ is the distance between the two communicating nodes, $l = \{A,G\}$ for A2G or G2G links, and $n = \{L,N\}$ for links in LoS or NLoS. The parameters $\beta_{ln}$ are the floating intercepts while the parameters $\alpha_{ln}$ are the path loss exponents. In order to avoid convergence issues with the numerical integration of the expressions for G2G links, as in \cite{Khoshkholgh2019}, the distance $r$ is replaced by $1+r$ in the path loss model.\smallskip

\paragraph{Small-scale fading}
We assume that the propagation is affected by Rician fading. Knowing the ground distance $r$ (in two dimensions) between the two communicating nodes, it is modelled by a random variable $h_{ln|r}$, with normalised power and a distance-dependent K factor $K_{ln}(r)$, different for A2G/G2G links in LoS/NLoS. We suppose that the small-scale fading of different links are independent, and that the random variables at different time instants for the same link are independent.\smallskip

\paragraph{LoS probability (G2G links)}
In this paper, the 3rd Generation Partnership Project (3GPP) $d_1/d_2$ model is selected, with a LoS probability given by \cite{Sun2015}
\begin{equation}
p_{GL}(r) = \min\left(\frac{d_1}{r},1\right) \left(1 - \exp{-\frac{r}{d_2}}\right) + \exp{-\frac{r}{d_2}}, \tag{7}
\end{equation}
where $r$ is the ground distance. The parameters $d_1$ and $d_2$ are respectively called the near-field and far-field critical distances.\smallskip

\paragraph{LoS probability (A2G links)} A commonly used LoS probability function for A2G links is given by \cite{Al-Hourani2014}
\begin{equation}
    p_{AL}(r) = \frac{1}{1 + a\  \exp{-b\left(\frac{180}{\pi}\arctan\left(\frac{H_R}{r}\right)-a\right)}},
    \label{eqn:LOSprobA2G} \tag{8}
\end{equation}
where $r$ is the ground distance, and $a$ and $b$ are two fitting parameters, directly linked to three other parameters with a physical meaning \cite{Al-Hourani2014}: the ratio of built-up land area to the total land area, the mean number of buildings per unit area, and a scale parameter that describes the buildings' heights distribution according to a Rayleigh Probability Density Function (PDF).\smallskip

\paragraph{Noise} We suppose that every link is affected by a white Gaussian noise of constant power $\sigma^2$.
\section{Introduction}
\input{Introduction.tex}

\section{System Model}
\input{System_Model.tex}

\section{Coverage Probability}
\label{section:coverage_probability}
\input{Coverage_Probability.tex}

\section{Numerical Results}
\input{Numerical_Results.tex}

\section{Conclusion}
\label{section:conclusion}
\input{Conclusion.tex}

\bibliographystyle{ieeetr}


\section{Introduction}
\label{sec:intro}
COVID-19 has wreaked havoc around the world and the spread of COVID-19 has not been curbed totally by June, 2020. According to the situation report of the WHO, there have been 9,653,048 confirmed cases and 491,128 deaths by 27 June 2020\cite{world2020coronavirus}. Recently, many researchers have developed a variety of detection methods of COVID-19, such as reverse transcription-polymerase chain reaction (RT-PCR) technology\cite{lan2020positive, hussein2020point}, X-ray\cite{wang2020covid} and chest CT imaging\cite{bernheim2020chest}. In addition to the need for breakthroughs in testing techniques, SARS-CoV-2 must be controlled at source to prevent a sustaining surge in COVID-19 confirmed and death cases. Small respiratory droplets can transmit SARS-CoV-2, which is known to be transmissible from presymptomatic and asymptomatic individuals\cite{howard2020face}. One thing to reduce virus's spreading is wearing masks in public\cite{howard2020face}. Advice on using face masks proposed by the WHO\cite{10665-332293}, analysis\cite{howard2020face} and researches\cite{greenhalgh2020face} have supported that wearing masks can prevent the spread of potential viruses during pandemic. An evidence research even thinks that wearing masks in public is a symbol of social solidarity among the response to the worldwide pandemic \cite{cheng2020wearing}. The WHO recommends that the surgical face masks used by health workers should be discarded once used and during severe shortages the extended use of faces masks can be considered \cite{world2020rational}. Because face masks only work when being worn within their periods of validity. On the other hand, more often in daily life, the majority of people wear surgical face masks in short-time but high-frequent way and in low-risk or middle-risk environment. The face masks may be used twice or three times, even four or five times if in relatively safe environment and maintained properly, which is reasonable and common especially under the condition where face masks are in shortage.

In this case, there is no doubt that the service life of masks depends on the use methods, use intensity, environmental condition and the material of masks, etc. In daily life, it often happens that common people forget the service stage of the current mask after having used it in the high-frequent but short-time way. Meanwhile, the service life of face masks decreases differently after being used in different conditions. As a result, the masks' remaining time of effective protection varies. With those faces masks whose use time are unknown, continuing to use them may lead to poor protection and other unexpected incidents; discarding them directly is an inappropriate behavior. It is a waste of face masks, especially in the case that masks are in shortage during the epidemic. It is necessary to detect the service stage of our face masks, but there are few researches on the service stage of face masks. Therefore, we aim at studying what service stage the face mask is in, which may play a positive role in protecting the uninfected people and reducing the spread of virus.

According to Javid et al., wearing a mask in public may become our unified action in the fight against COVID-19\cite{javid2020covid}. Due to the short service life and worldwide use of masks, it is not practical for masks to adopt special machines to centralized detection. Hence, it is urgently needed to develop an easy-to-operate, portable testing device for detecting anytime and anywhere. Therefore, in this research, we propose a portable mask service stage detection system based on mobile phone. Its work procedures are as follows: 1) users first take a micro-photo of the mask being used with a mobile microscope; 2) then the micro-graph of the mask is uploaded in the WeChat Applet; 3) the result of back-end detection eventually can be obtained through WeChat Applet. This detection device is simple, easy to operate and effective.

In the common fight, the ordinary people often wear surgical masks while N95 masks are only recommended for healthcare workers or professionals at high risk of coming into contact with patients\cite{long12381effectiveness}. Radonovic et al. analyzed the data which came from 7 health care delivery systems and 4 seasons of peak viral respiratory illness, and found that in this trial, wearing N95 respirators and wearing medical masks resulted in no significant difference in the incidence of laboratory-confirmed influenza\cite{radonovich2019n95}. Additionally, the research suggests that N95 respirators should not be recommended for the general public, namely, those not in close contact with influenza or suspected patients\cite{long12381effectiveness}. Using N95 respirators may lead to discomfort, like headaches\cite{cowling2010face}. Even a previous study \cite{chen2017herd} has reported that there is an inverse relationship between the level of wearing N95 respirators and the risk of clinical respiratory illness. Therefore, we choose surgical mask as the research object, which is the most widely used. It makes experimental data more scientific and ultimately our detection system can benefit more people to the greatest extent. In this study, after obtaining the micro-graphs of masks at different wearing periods, the texture features of the photos are analyzed.

Methods of texture analysis include fractal analysis, Fourier transformation and the gray level co-occurrence matrix (GLCM)\cite{kanai2020discriminant}. Haralick et al. reported that the GLCM was a method to quantify the spatial relationship between adjacent pixels in an image\cite{haralick1973textural}. GLCM is widely used in disease detection \cite{gibbs2003textural,rathore2014ensemble}, skin texture analysis\cite{ou2014vivo}, and defect detection\cite{shabir2019tyre}, etc. Given the above, GLCM is feasible and efficient in image texture feature analysis. Therefore, we choose GLCM to analyze the texture features of masks.

Due to K Nearest Neighbor (KNN) method's simple implementation and distinguished performance in classification tasks\cite{zhang2017learning}, it is widely employed in text categorization\cite{du2019parallel,li2019improved}, medical domain\cite{hamed2020accurate}, eating patterns exploration\cite{newby2004empirically}, leaf disease identification\cite{krithika2017individual}, indoor localization\cite{hoang2018soft}, etc. Hamed et al. compared the performance of KNN variant (KNNV) algorithm (which derived from KNN) and other three algorithms in the classification of COVID-19 patients, and found that KNNV algorithm could classify COVID-19 patients more efficiently and accurately\cite{hamed2020accurate}. On the basis of the texture features extracted from the micro-photos taken in different period time, this paper employs KNN algorithm to establish the classfiers for detecting the use time of mask.

In this paper, we propose a portable face mask service stage detection system, which is easy to operate and can work anytime and anywhere. It is based on the texture features of face masks and ground-truth of mask service stage data. To further quantify the differences of masks in different periods, the method of texture analysis is introduced. After extracting texture features from the micro-photos of face masks in different use periods, KNN algorithm is performed to detect which service period the face mask belongs to. The testing data are afterward used to validate the effectiveness of the proposed system.

The main contributions of this paper are threefold. First, we combine the texture features of face masks' micro-photos and their corresponding service stages based on the method of texture analysis. Based on GLCM, the texture features are successfully extracted from micro-photos obtained by high magnification lens attached to the mobile phone. Subsequently, we propose a detection method to judge the service stage of face masks with a KNN algorithm. Finally, based on the two contributions mentioned above, we have implemented a detection Wechat Applet for face masks using the data collected from three persons' wearing face masks in daily life, which may contribute to protecting the uninfected people and reducing the spread of virus.

\section{Methods}
\label{sec:method}

A brief introduction to the proposed face mask service stage detection method is shown below. We initially use the portable and delicate mobile microscope device to get the micro-photos of face masks. After obtaining the picture, the first step is to extract texture features data from it. During the extraction process, we use GLCM to capture the contrast, energy, correlation, and homogeneity of the object face mask. Subsequently, KNN algorithm is employed to propose a classification model with the input texture features and the corresponding ground-truth use time. Eventually, we package the system into a library function in Matlab and embed it to the Wechat Applet we develop. This has been proved effective by similar previous researches\cite{jiang2020detection,imran2020ai4covid}.

\subsection{Overview of Data Acquisition}
\begin{figure}[ht]
	\centering
	\includegraphics[width=0.45\textwidth]{F01.png}
	\caption{The overview of the portable detection device's ``Hardware''. The ``Hardware'' is composed of a mobile microscope, a phone with available system camera, etc. Specifically, a) and b) represent the mobile microscope and system camera, respectively; and c) shows the final removable device mainly assembled by a) and b).}\label{F01}
\end{figure}

\begin{figure}[ht]
	\centering
	\includegraphics[width=0.5\textwidth]{F02.png}
	\caption{Scenario of the portable detection device's working. While the device is on work, the flash must meet the condition that the flash is always ``on'' to provide enough light.}\label{F02}
\end{figure}

\begin{figure}[ht]
	\centering
	\includegraphics[width=0.5\textwidth]{F13.png}
	\caption{The three inside photographed locations of the face mask: left location (about one-third of the face mask and close to the left), middle location (about two-thirds of face mask), right (about one-third of the face mask and close to the right).}\label{F13}
\end{figure}

Fig. \ref{F01} displays the hardware configuration of the detection device: a mobile microphone (TIPSCOPE, Convergence(Wuhan) Technology Ltd., China), a phone with available system camera, etc. The detachable mobile microphone, which is the key to the hardware, is particularly chosen to display the details of face masks viz. ash particle, droplets or other debris more explicitly. It should be noted that the maximum magnification can reach 400x if the system camera and mobile microscope work together. In this research, we use system camera without magnification to get more scope of face masks.
Fig. \ref{F02} shows how to collect the experimental data. To ensure the validity of the collected data, we photograph three inside locations of the face mask each time: left location (about one-third of the face mask and close to the left), middle location (about two-thirds of face mask), right (about one-third of the face mask and close to the right), as shown in Fig. \ref{F13}. In order to contain more usage scenarios, the face masks are worn in different conditions. Meanwhile, the flash of the phone's camera should be always ``on'' to provide enough light. We collect the images of a face mask used from day zero (new) to day five.
Based on the conditions mentioned above, three images are immediately photographed after the face mask has been used for an arbitrary hour of the daytime.

\subsection{Gray Level Co-Occurrence Matrix (GLCM)}
Texture can be used to characterize the tonal or gray-level variations in an image\cite{hall2017glcm}. The gray level co-occurrence matrix (GLCM) is chosen to disclose the texture features hiding in the face masks. GLCM is a second-order statistical method which calculates the frequency of pixel pairs with the same gray-level in an image and uses additional knowledge obtained from spatial pixel relations\cite{wu2015active}. The co-occurrence matrix uses edge information to embed the distribution of gray-scale transformation\cite{xing2019multilevel}. Owing to the fact that much of the information required is embedded in GLCM, it emerges as a simple but effective technique.

14 measures of textural features are proposed by Haralick et al., and some of these measures are related to specific textural characteristics of the image\cite{haralick1973textural}. In our study, we choose contrast, correlation, energy and homogeneity as measures. Specifically, contrast refers to the drastic change in gray level between adjacent pixels. High contrast images own high spatial frequencies. Correlation represents the linear correlation in an image. The higher the correlation value, the more linear the gray-scale relationship between adjacent pixel pairs. Energy stands for texture uniformity or pixel pair repetitions. High energy is produced when the distribution of gray level values is constant or periodic. Homogeneity is sensitive to the existence of near diagonal elements in a GLCM, indicating the similarity of adjacent pixels in gray scale.

With the increase of the use times, the wear of the mask causes a change in micro-structure of the mask's some parts, which leads to the difference between the worn-out parts and the unworn parts. Thus, the images in different the mask's service stage
may perform variously in correlation and energy. Meanwhile, due to humans' respiration, cough and other respiratory related activities, more and more small particles impurities and droplets will adhere to the inside surface to the face mask. They probably result in the difference in the four measures mentioned above, among the images of different service stages.

Every element in GLCM contains second-order statistics, probability values for changes between gray levels $i$ and $j$ for a particular displacement $d$ and angle $\theta$, labeled as $p(i,j)$ (normalized)\cite{lofstedt2019gray}. Let $M$ be the number of gray levels in the image, thus, the size of GLCM is $M\times M$. In this research, we set $M$ as $8$. Subsequently, the measures are calculated as

\begin{equation}\label{E01}
Contrast=\sum_{i=0}^{M-1} \sum_{j=0}^{M-1} (i-j)^2 \cdot p^2(i, j)
\end{equation}

\begin{equation}\label{E02}
Homogeneity=\sum_{i=0}^{M-1} \sum_{j=0}^{M-1} \frac{1}{1+(i-j)^2} \cdot p(i, j)
\end{equation}

\begin{equation}\label{E03}
Energy=\sqrt{\sum_{i=0}^{M-1} \sum_{j=0}^{M-1} p^2(i, j)}
\end{equation}

\begin{equation}\label{E04}
Correlation=\sum_{i=0}^{M-1} \sum_{j=0}^{M-1} (i-\mu)\cdot (j-\mu)\cdot p(i, j)/\sigma^2
\end{equation}

where $\mu$ is the mean of GLCM while $\sigma^2$ stands for the variance of GLCM.

\subsection{K Nearest Neighbor (KNN) Algorithm}
Considering that this task is relatively simple, so the model can be constructed with a relatively simple pattern recognition method. At the same time, in view of the deployment to the mobile phone terminal, the established model should be easy to deploy to the mobile phone terminal and convenient to maintain later. Compared with other algorithms, due to its simple implementation and significant classification performance, KNN is a very popular method in statistics and ranks top ten data mining algorithms{\cite{zhang2018efficient, 7401007, wu2007top,zhang2012nearest}.

KNN algorithm is different from model-based methods which first use training samples to build a model, and then predict the test samples through the learned model\cite{yu2014a, zhu2015targeting,shao2014feature}. No training phase is required for the model-free KNN method. Instead, it conducts classification tasks obeying the following procedures: first, training samples are attached labels; next, the distance is calculated between the test sample and the training samples in each label; subsequently, after comparing the distances and obtaining the test's nearest neighbors, the serial number is obtained; finally, classification is generated, as shown in Fig. \ref{F04}. The distance is computed as:
\begin{equation}\label{E05}
	Distance=\sqrt{\sum_{j=1}^{n} (test\ sample-training\ sample_{j})^2}
\end{equation}
where the $n$ stands for the number of labels.

\begin{figure*}[ht]
		\centering
		\includegraphics[width=0.9\textwidth]{F04.png}
		\caption{The flowchart of KNN algorithm. First, labels are attached to the training  samples. For every test sample, the $distance_{i}$ is calculated between the test sample and the training samples in each label, which follows the equation \ref{E05}. By comparing the distances, the test sample's nearest neighbors viz.the minimum distance are obtained, and serial number $i_{min}$ is then produced. Eventually, the final classification is set as label $i_{min}$.}\label{F04}
\end{figure*}

The main task in our research is to train the optimal $k$ value after carrying out KNN classification successfully according to the efficiency of the classification performance.

\subsection{Portable device based on mobile phone for detecting use time of mask}
After optimizing the system, the system is embedded to the Wechat Applet we have developed. To enrich our experimental data, we enlarge the data collection time from the first day to the fifth day (about an hour used time each day). It is generally recommended to replace the surgical mask every four hours. If the face mask is used for a long time, the large particles will be blocked on the mask surface or the ultrafine particles will be blocked in the pores of the mask filter material, resulting in the decrease of filtration efficiency and the increase of respiratory resistance\cite{shiyong}. According to the detection time, the detection results are within three categories viz. type I: normal use (day 0 to 1 viz. the face masks can be used securely), type II: early warning (day 2 to 3 viz. the face masks can be used safely, but it is close to the end of its service life) and type III: not recommended (day 4 to 5 viz. the face mask can be used in shortage, but it is better to change a new one). We choose the `day' as unit based on the fact that we use face masks in high-frequent but short-time way and ordinary people hardly use face masks for 4 hours continuously. Additionally, on basis of using the same face mask several times, we classify the face masks belonging to day 4 and day 5 the type `not recommended' to counterbalance the increase in use times.
\begin{figure*}[ht]
	\centering
	\includegraphics[width=1\textwidth]{F14.png}
	\caption{The workflow chart of the face mask detection system. a) stands for collecting the micro-image of the face mask being used. b) with its branches b1), b2) and b3) means registering an account used to receive the detection result and the following uploading task. c) is the result sent by the Wechat Applet from the back-end. (The face mask used in the sample is in its service life of type I: normal use.)}\label{F03}
\end{figure*}
Fig. \ref{F03} reveals the workflow of the detection. After we have obtained the micro-photo of the face mask with mobile microscope and registered an account in the Wechat Applet, the micro-image is uploaded to the back-end of the system. Accordingly, the detection result will be sent in the message from the back-end.

\subsection{Evaluation Metrics}
To verify the efficiency of the detection system, macro-measures viz. macro precision, macro recall, and macro $F_{1}$ are considered.

1) Confusion matrix: we assume that ``Positive'' means the positive samples and ``Negative'' means the negative samples. Meanwhile, ``True'' represents that the prediction is right while ``False'' represents that the prediction is wrong. As a result, ``TP'' and ``TN'' mean that the positive sample is classified as ``Positive'' and the negative sample  is labeled as ``Negative'', respectively. ``FP'' and ``FN'' represent that the negative sample is labeled as ``Positive'' and the positive sample is classified as  ``False''. The four indicators make up the confusion matrix.

2) Precision: precision is only used to evaluated the classification ability of the positive samples within the range from 0 to 1. It is obvious that the larger precision is, the more effective the system is. It is computed by:
\begin{equation}\label{E07}
Precision=\frac{TP}{TP+FP}
\end{equation}

3) Recall: it is a ratio from 0 to 1. Obviously, the more it is close to 1, the better the system is. The calculation equation is:
\begin{equation}\label{E08}
Recall=\frac{TP}{TP+FN}
\end{equation}

4) $F_{1}$: it is a harmonic mean of recall and precision. In this study, we consider the weight of recall and precision the same, which means attaching the weight of 0.5 to either of them. It is calculated by:
\begin{equation}\label{E09}
F_{1}=\frac{2*Precision*Recall}{Precision+ Recall}
\end{equation}

7) Macro-measures: metrics are calculated for each label, and then the weighted means are produced (normally the weight is the same except that the number of the samples in each label varies greatly).
\begin{equation}\label{E13}
Macro\_Precision=\frac{1}{n}\sum_{i=1}^{n}Precision_{i}
\end{equation}

\begin{equation}\label{E14}
Macro\_Recall=\frac{1}{n}\sum_{i=1}^{n}Recall_{i}
\end{equation}

\begin{equation}\label{E15}
Macro\_F_1=\frac{1}{n}\sum_{i=1}^{n}Macro\_{F_1}_{i}
\end{equation}
where $n$ means the total number of labels, and the $Macro\_{F_1}_{i}$ is $F_{1}$ of every label calculated by Equation \ref{E09}.
\section{Experimental Results and Analysis}

In this study, we have collected 87 micro-images (excluding the invalid data) of surgical masks via the mobile microscope. The life span of face masks have little relationship with different people but much with the condition of use environment. Meanwhile, people's doing when wearing face masks also influences the life span of face masks. As a result, the data cover the conditions: speaking frequently, talking with others, running, shopping in supermarkets, shopping in vegetable markets, riding and wandering, etc.

Every time three micro-photos (the left, the middle and the right) are photographed, they are initially averaged to eliminate existing disturbance, and then one final training sample is produced. As a result, after extracting texture features, 29 sets of training samples are obtained. The original 87 micro-images are used to carry out validation experiments. The evaluation indicators mentioned above are elaborated as the follows.

\subsection{Experimental Results}
\begin{figure}[ht]
	\centering
	\includegraphics[width=0.45\textwidth]{F05.png}
	\caption{The confusion matrix of the classified results. Each row is the number of real labels and each column is the number of classified labels. Type I, Type II and Type III represent `normal use', `early warning' and `not recommended', respectively.}\label{F05}
\end{figure}
To figure out the detailed information about the classification of the use time, we plot the confusion matrix of the three types as demonstrated in Fig. \ref{F05}. As can be seen from the results, 71 samples among 87 samples are predicted correctly, which suggests the model performs relatively well on the whole.

\begin{table}[htbp]
	\centering
	\caption{The Macro Measures of Type I, II and III.}
	\begin{tabular}{rrcccc}
		\toprule
		&       & \multicolumn{2}{c}{Macro Precision} & \multicolumn{2}{c}{82.87\%$\pm$8.50\%} \\
		\multicolumn{2}{c}{Macro measures} & \multicolumn{2}{c}{Macro Recall} & \multicolumn{2}{c}{81.98\%$\pm$7.50\%} \\
		&       & \multicolumn{2}{c}{Macro $F_1$} & \multicolumn{2}{c}{81.66\%$\pm$1.40\%} \\
		\bottomrule
	\end{tabular}%
	\label{tab1}%
\end{table}%
According to the equations of macro-measures, we can easily see that macro-measures are convincing and scientific. The macro-precision of the system reaches 82.87\%$\pm$8.50\%, which proves the efficiency of the detection system.
\subsection{Influence of K Factor on the Efficiency of System}
\begin{figure}[ht]
	\centering
	\includegraphics[width=0.45\textwidth]{F06.png}
	\caption{The influence of K factor on the performance of system. The Y axis on the right stands for the number of correctly predicted samples, and the Y axis on the left shows the corresponding accuracy. Intuitively, it can be seen that the optimal K is 6.}\label{F06}
\end{figure}
As mentioned above, we specially aim at obtaining the optimal $k$ to get the most effective detection system. In KNN algorithm, K plays a key role in the performance of system. In other words, if the K factor is too small, there may exist the following problems: 1) the phenomenon of over fitting; 2) being easy to influenced by outliers; and 3) the system is too complex. While the K factor is  too large, there may appear the elaborated problems: 1) the phenomenon of under fitting; 2) the system is constrained by the sample equilibrium; and 3) the system is too simple.

Hence, in this study, we analyze the correctly classified samples' number and accuracy under the condition that K's values vary from 4 to 16, as demonstrated in Fig. \ref{F06}. Accordingly, we can find the optimal-k-value is 6, with the the number TP 71 and the precision 82.87\%. Additionally, K going from 4 to 6, the curve is roughly go up. For K ranging from 6 to 16, there are two stages in the descent tendency of the curve. First, the curve declines rapidly when K is from 6 to 9. Then it oscillates downward from 9 to 15, and the performance of system descends sharply when the K is changed to 16. The phenomenon referred above conforms to the theoretic prediction. Surprisingly, it is exactly consistent with the conclusion that the fixed optimal-k-value for all test samples should be $k=\sqrt{n}$ (where $n>100$ and $n$ is the number of training samples proposed by Lall and Sharma\cite{95WR02966}.
\subsection{Comparison of measures in three types categories}

\begin{figure*}[ht]
	\centering
	\includegraphics[width=1\textwidth]{F07.png}
	\caption{The micro-photos of three types face masks and their corresponding GLCMs' measures. a), b) and c) represent the face masks in `normal use' stage, `early warning' stage and `not recommended' stage, respectively. d), e), f) and g) are curve of four measures viz. contrast, correlation, energy and homogeneity.} \label{F07}
\end{figure*}
we analyze the measures of all samples to inquire the tendency of modeling measures over time changes. Fig. \ref{F07} displays the micro-photos and their corresponding GLCMs' measures of three different service stages viz. `normal use' stage, `early warning' stage and `not recommended' stage. As we can see from the micro-photos, when the use time increases, the mask becomes more and more dirty and there are many impurities and small droplets in the holes of the mask. It proves the accuracy of the analysis mentioned above that the protection becomes less effective over time. The curves reveal the variation tendency of the measures. Specifically, contrast and correlation increase while energy and homogeneity decrease over time. It should be noted that this rule can't be applied to all data, since the detection result is determined by composite force of the four measures. As a result, there exists the situation that some data don't obey the rule.

\subsection{Comparison of Experimental Results using original photos and micro-photos, respectively}
\begin{figure}[ht]
	\centering
	\includegraphics[width=0.45\textwidth]{F08.png}
	\caption{The original photo and the micro-photo of the same face mask. a) displays the micro-photo obtained by microscope while b) shows the original photo gained only by mobile phone's camera.}\label{F08}
\end{figure}

To validate the necessity of microscope, we record the results of using the original photo and micro-photo, which are shown in Fig. \ref{F08}. The two photos represent the same location of the face mask and belong to the type `not recommended'. In Fig. \ref{F08}, we can see many details, such as droplets in the holes, etc. Meanwhile, we can see no details in the original photo except the obvious stain. It can be inferred that we will tell little difference if there is no stain or other obvious dirt. The detection result of the micro-photo is `not recommended' while the original photo 's result is  `normal use', which verifies the indispensability of the microscope in this system. In other words, the photos obtained without using microscope aren't able to be used as the detection input.
\subsection{Dichotomous Model of Experimental Results}

On the basis of system's working normally, we further have the idea whether we can simplify the system. As a result, we propose a dichotomous model with only two types viz. `normal use' and `not recommended'. Considering that `normal use' and  `not recommended' type face masks both can protect people effectively, we merge them into one type named `normal use'. Accordingly, the other type is `not recommended'. In the test samples, the number of `normal use' is 60 (54 samples are detected correctly) and `not recommended' is 27 (23 samples are detected correctly). To counterbalance the uneven distribution of samples, we assign the weight viz. $\frac{60}{87}$ and $\frac{27}{87}$ to each type. As expected, the accuracy becomes better and increases to 88.51\%. Since there are two types, it is easier for the model to tell which type the detected face mask belongs to. In some simple cases, we may consider this scheme.
\subsection{Limitation of face mask condition }
Besides the system mentioned above being able to be used to detect its service life, some conditions of the face mask itself may also affect their use life. For example, there may be lint inside the face mask when using face mask a few times. In that case, we can choose to discard the mask and use a new one to ease the discomfort. Meanwhile, if the face mask being used gets wet by accident, we must stop using it whatever stage it belongs to.

\section{CONCLUSION}
\label{sec3}

In this paper, we propose a service stage detection method based on a mobile microscope which can be used to obtain the micro-photos of the face mask being used.

In our detection method, we first get the micro-photos of face mask being used, which may reflect some details such as droplets or other obvious dirt. Then, we extract texture features from the micro-photos using the method of GLCM and choosing four measures viz. contrast, correlation, energy and homogeneity as the features. Subsequently, KNN method is applied to work on the detect the service stage. In validation experiments, the obtained system achieves a relatively good result with a precision of 82.87\%$\pm$8.50\% on the testing dataset. The result of our experiments provide an idea of using photos to detect the use time of face masks to distinguish whether the face mask can be used sequentially or not. Faced with the current severe situation of COVID-19 and possible shortage of face masks, our research may work as an assistant to help the common people to use face masks more correctly. Furthermore, it can exert positive influence on protecting uninfected people and stopping the possibly infected patients from spreading virus.

In future research, based on the current portability and practicability, we will try using a more stable algorithm to reach the goal of achieving a higher accuracy. Besides, different types of face masks should be taken into consideration to expand the types of detecting objects.

\small
\bibliographystyle{IEEEbib}


\section{Introduction}

The core of Neutron Stars (NSs) reaches densities beyond that of nuclear saturation, allowing a phase transition from hadronic to Strange Quark Matter (SQM), which is speculated to be the truly ground state of strong-interacting matter (Bodmer-Witten's conjecture)~\cite{Bodmer, Witten}. Experimentally, the quark--gluon plasma, in the limit of high temperatures and low densities, has been explored in the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory (BNL) and in the Large Hadron Collider (LHC) at CERN. Meanwhile, the opposite regime of low temperatures and high densities will be studied in two experimental facilities that are currently under construction, the Facility for Antiproton and Ion Research (FAIR) at GSI and the Nuclotron-based Ion Collider Facility (NICA) at JINR. Nevertheless, since NSs cores might be the natural habitat for SQM in the latter regime, the properties of matter under those conditions could also be inferred through astronomical observations.

A subclass of Neutron Stars composed by a gas of degenerated quarks was first proposed by Itoh in $1970$ \cite{Itoh}, and since then many microscopic models have been proposed for those objects.  The main difference between models are how the strong interaction is described and  what type of SQM phases is considered inside the star \cite{Camenzind}.
Compact objects with quark matter could then be either SSs composed entirely by SQM, or hybrid stars when there is a core of quark matter surrounded by a layer of hadrons.
SSs models are interesting because they describe compact objects with maximum masses around $1.5 \, M_{\odot}$ and radius of $4-8$~km that could account for some observations of cold, dense and small Compact Objects (COs) that do not fit the standard NSs models (see Table \ref{candidatas}) \cite{Lattimer2003, Abubekerov, Taparati, Guver2010_1, Guver2010_2, Elebert, Rawls}. On the other hand, Hybrid Stars have received a lot of attention  because they can reach $2 M_{\odot}$, which is a robust observational constrain on NSs mass \cite{Demorest}.  Recently, these models have gained a new support due to the works by Eemeli Annala et al.\cite{Annala}, where a wide set of theoretical EoS from particle and nuclear physics is compared with benchmark results stemming from Gravitational Waves (GWs) measurements of NSs collisions. According to their results, a star with mass $\sim 2 \, M_{\odot}$ and radius $\sim 12$~km is more likely to have a quark core of approximately $6.5$~km, than to be formed exclusively by baryons.
\begin{table*}[ht!]
	\caption{Masses ($M$) and  radii ($R$) for observed candidates of SSs.}
   	\label{candidatas}
	\centering
	\begin{tabular}{|l|l|l|}
		\hline
		Observed Object & $M(M_{\odot})$  & R(km)\\

		\hline
		\hline
		RXJ185635-3754 \cite{Lattimer2003} & $1.20\pm0.2$ & $8.00\pm1$  \\
   		HerX-1 \cite{Abubekerov, Taparati} & $0.85\pm0.15$ & $8.1\pm0.4$ \\
   		4U1608-52 \cite{Guver2010_1} & $1.74\pm0.14$ & $9.3\pm1$ \\
   		4U1820-30 \cite{Guver2010_2} & $1.58\pm0.06$ & $9.1\pm0.4$ \\
   		SAXJ1808.4-3658 \cite{Elebert} & $0.9\pm0.3$ & $7.9\pm1$ \\
   		Vela X-1 \cite{Rawls} & $1.77\pm0.08$ & $9.56\pm0.08$ \\
   		4U1538-52 \cite{Rawls} & $0.87\pm0.07$ & $7.87\pm0.21$ \\
   		SMC X-1 \cite{Rawls} & $1.04\pm0.09$ & $8.3\pm0.2$ \\
   		LMC X-4 \cite{Rawls} & $1.29\pm0.05$ & $8.83\pm0.09$ \\
   		Cen X-3 \cite{Rawls} & $1.49\pm0.08$ & $9.18\pm0.13$ \\
   \hline
	\end{tabular}
\end{table*}

The study of SSs starts by seeking the SQM Equation of State (EoS). Inside SSs quarks might be unpaired or paired forming a color superconductor state \cite{Ferrer}. However, so far there are neither \textit{ab initio} nor perturbative QCD calculations that leads to the desired EoS in the regime of high densities and zero temperature\cite{Buballa:2003qv}. Consequently, SQM inside COs is usually described through phenomenological models that mimic the main features of QCD \cite{Buballa:2003qv}. One of the most used is the MIT--Bag model \cite{Chodos}, where quarks are considered as quasi-free particles confined into a \textquotedblleft bag\textquotedblright and having fixed masses. This model reproduces confinement and asymptotic freedom with the use of only one external parameter, the bag energy $\mathtt{B_{Bag}}$.

On the other hand, Compact Objects have extreme magnetic fields \cite{Shapiro}. The observed surface magnetic fields in Neutron Stars are in the range of $10^{9}-10^{15}$ G \cite{1993ApJ...408..194T, 1998Natur.393..235K}, while their inner magnetic fields are estimated to be as high as $5\times 10^{18}$ G \cite{Daryel2015, Lattimerprognosis}. As it is well known, the energy momentum tensor of matter under a magnetic field is anisotropic \cite{Ferrer}, a feature that can be interpreted as the gas exerting different pressures in the parallel and perpendicular directions with respect to the magnetic axis \cite{Bali:2014kia}, and that leads to non-spherical stars \cite{Konno:1999zv,Chatterjee:2014qsa,Rizaldy,prd}.
To model the structure of magnetized COs, we have previously used Tolman-Oppenhimer-Volkoff (TOV) equations independently for each pressure \cite{Felipe2008, Felipe2009JPhG,ARD,Perez2019}, and an approximate solution of Einstein's equations for an axially symmetric metric in cylindrical coordinates\cite{Daryel2015}. This leads to the description of two different stellar sequences (one for each pressure), which prevents  to calculate the total mass of the star. That is why, as an attempt to properly describe the macroscopic structure of magnetized COs,  we derived a set of TOV-like structure equations from an axially symmetric metric in spherical coordinates, the $\gamma$-equations \cite{prd}, that allow to describe spheroidal objects as long as their shape is nearly spherical.  These equations  have been used to model White Dwarfs \cite{prd}, the hypothetical magnetized BEC stars~\cite{Gretel} and, preliminary, the Strange Stars \cite{Perez2019}.

In this paper we returned to SSs \cite{Perez2019,Felipe2009JPhG,Felipe2008}, emphasizing in the magnetic field effects on the stability of the SQM and the spheroidal stellar configurations. In addition, we compute the eccentricity, the redshift and the mass quadrupole moment and compare them with those of spherical Strange Stars. This study is essential and is also the starting point of a more ambitious project related to the study of magnetized Hybrid Stars.

In Section \ref{sec2} we analyze the stability of magnetized SQM in astrophysical conditions and its dependence on the density, the bag energy and the magnetic field. Besides, we revisited the magnetized  EoS for SSs and discuss the magnetic field effects on the energy density and pressures. The $\gamma$-equations are presented in Section \ref{sec3} with their corresponding numerical results for magnetized Strange Stars. Section \ref{sec4} is devoted to the computation of the redshift and the mass quadrupole moment. Concluding remarks are given in Section \ref{sec5}.

\section{EoS of magnetized Strange Stars}\label{sec2}

We consider a Strange Star composed by SQM and electrons under the action of a uniform and constant magnetic field oriented in the $z$ direction ${\bf B} = (0,0,B)$. As pointed out before, we use the
 phenomenological MIT--Bag model \cite{Chodos}, where quarks are assumed as quasi-free particles confined into a \textquotedblleft bag\textquotedblright\, that reproduces the asymptotic freedom and confinement through the B$_{\text{bag}}$ parameter, a binding energy which is added to the quarks energy and subtracted from their pressure \cite{Chodos}.

For a magnetized gas of quarks and electrons, the pressure and the energy density are obtained from the thermodynamical potential \cite{Felipe2008}
\begin{equation}\label{Thermo-Potential}
\Omega_f(B,\mu_f,T)  = -e_fd_f B T\int_{-\infty}^{\infty} \! \frac{dp_3}{4\pi^2}  \sum_{l=0}^{\infty}g_l
  \times\sum_{p_4}\ln\bigg[(p_4+i\mu_f)^2+ \varepsilon^2_{lf} \bigg],
\end{equation}
where $l$ stands for the Landau levels and $f=e, u, d, s$ for the electrons and each quark flavour; $d_f$ is the flavour degeneracy factor\footnote{$d_e=1$ and $d_{u,d,s}=3$.} and $g_l=2-\delta_{l0}$ includes the spin degeneracy of the fermions for $l\neq0$. Moreover, $T$ is the absolute temperature, while $\mu_f$, $m_f$ and $e_f$  are the chemical potential, the mass and the charge of each particle respectively. The spectrum of charged fermions coupled to a magnetic field is $\varepsilon_{lf}=\sqrt{p_3^2+2|e_fB|l+m_f^2}$ \cite{Felipe2008}.

Eq.~\eqref{Thermo-Potential} can be divided in two contributions
\begin{equation}\label{Thermo-Potential-2}
\Omega_f (B,\mu_f,T)=\Omega^\text{vac}_f(B,0,0)+\Omega^\text{st}_f(B,\mu_f,T),
\end{equation}
with
\begin{eqnarray}
\Omega^\text{vac}_f(B,0,0) &=&-\frac{e_fd_f B}{2\pi^2}\int_{0}^{\infty}\hspace{-2mm}dp_3\sum_{l=0}^{\infty}g_l\varepsilon_{lf} \label{vaccon},\\\label{omegast}
\Omega^\text{st}_f(B,\mu_f,T) &=& -\frac{e_fd_f BT}{2\pi^2}\int_{0}^{\infty}\hspace{-2mm}dp_3\sum_{l=0}^{\infty}g_l
\ln\hspace{-1mm}\big[1+e^{-\beta(\varepsilon_{lf}\pm\mu_f)}\big]
\label{stcon}\hspace{-0.5mm}.
\end{eqnarray}
 The vacuum contribution Eq.~\eqref{vaccon} does not depend on the chemical potential nor on the temperature and presents an ultraviolet divergence that must be renormalized \cite{Schwinger}, resulting in
\begin{equation}\label{vaccon2}
\Omega^\text{vac}_f(B,0,0)= \frac{d_fm_f^4}{24\pi^2}\left(\frac{B}{B^c_f}\right)^2\text{ln}\frac{B}{B^c_f}.
\end{equation}
In Eq.\eqref{vaccon2} $B^c_f=m^2_f/e_f$ is the  critical magnetic field\footnote{The magnetic field at which the cyclotron energy of the particles is comparable to their rest mass. For electrons this is the so-called Schwinger field.}. For electrons $B^c_e\sim10^{13}$~G while for quarks up, down and strange we have $B^c_u\sim10^{15}$ G, $B^c_d\sim10^{16}$ G and $B^c_s\sim10^{19}$ G, respectively.

COs have temperatures much smaller than the Fermi temperature of the gases that compose them \cite{Camenzind}. Hence, a good approximation is to compute the thermodynamical potential of these gases in the degenerate limit (T$\rightarrow0$) \cite{Shapiro, Camenzind}. In that case, the statistical part of the thermodynamical potential becomes
\begin{equation}\label{stcon2}
\Omega^\text{st}_f(B,\mu_f,0)=-\frac{e_fd_f B}{2\pi^2}\int_{0}^{\infty}dp_3\sum_{l=0}^{\infty}g_l\Theta(\mu_f-\varepsilon_{lf}),
\end{equation}
where $\Theta(\zeta)$ is the unit step function. From  Eq.~\eqref{stcon2}, we obtain
\begin{equation}\label{stcon3}
  \Omega^\text{st}_f(B,\mu_f,0)=-\frac{e_fd_f B}{4\pi^2}\sum_{0}^{l_{max}}g_l\bigg[\mu_fp_f^l
  -(\varepsilon_f^l)^2\text{ln}\bigg(\frac{\mu_f+p_f^l}{\varepsilon_f^l}\bigg)\bigg],
\end{equation}
where $p_f^l=\sqrt{\mu_f^2-(\varepsilon_f^l)^2}$ is the fermi momentum of the particles, $\varepsilon_f^l=\sqrt{m_f^2+2qBl}$ their ground state energies and $l_{\text{max}}=I[\frac{\mu_f^2-m_f^2}{2e_fB}]$ the maximum number of occupied Landau levels for fixed magnetic field and chemical potentials. $I[z]$ denotes the integer part of $z$.

Due to the high fermionic densities in the stars interior, the vacuum contribution Eq.~\eqref{vaccon2} is negligible with respect to the statistical one,  Eq.~\eqref{stcon3}~\cite{Ferrer}.
Therefore, the thermodynamical potential of the degenerate fermion system can be approximated to $\Omega_f(B,\mu_f,0)=\Omega^\text{st}_f(B,\mu_f,0)$.

Strange Quark Matter inside the star must be in stellar equilibrium \cite{Felipe2008}. So, we impose $\beta$ equilibrium, charge neutrality and baryon number conservation  to the system in terms of the particle densities $N_f=-\partial\Omega_f/\partial\mu_f$ and the chemical potentials. These conditions read \cite{Felipe2008}
\begin{subequations}\label{eqsqm}
    \begin{eqnarray}
    \mu_u+\mu_e-\mu_d=0 \,\,\,, \,\,\, \mu_d-\mu_s&=& 0,\\
    \sum_f  e_fN_f=0 \,\,\,, \,\,\

    \sum_{i=u,d,s}N_i=3N_B.
    \end{eqnarray}
\end{subequations}

With these considerations, the magnetized SSs EoS are
\begin{subequations}\label{EoSSS}
    \begin{align}
    E &= \sum_{f}\left[\Omega_f+ \mu_fN_f\right]+B_{\text{bag}} +\frac{B^2}{8\pi},\label{EoS1}\\
    P_{\parallel} &=- \sum_{f}\Omega_f-B_{\text{bag}}-\frac{B^2}{8\pi}, \label{EoS2}\\
    P_{\perp} &= -\sum_{f}\left[\Omega_f+B\mathcal{M}_f\right]-B_{\text{bag}}+\frac{B^2}{8\pi},\label{EoS3}
    \end{align}
\end{subequations}
where $\mathcal{M}_f=-\partial\Omega_f/\partial B$ is the magnetization. In Eqs.~\eqref{EoSSS} we can easily identify three different contributions. The first one is given by the sum over the thermodynamical quantities of the species and corresponds to the statistical contribution of each kind of particle. In the case of the perpendicular pressure Eq.~\eqref{EoS3}, this term includes a contribution that comes from the particle magnetization: $-B\mathcal{M}_f$ \cite{Chaichian}. The second terms in the EoS, $\pm B_{\text{bag}}$, are the ones ensuring asymptotic freedom and confinement for quarks \cite{Witten,Buballa:2003qv}. Finally, the last terms in Eqs.~\eqref{EoSSS} are the magnetic field pressures and energy density $P^B_{\perp}=E^B=-P^B_{\parallel}=B^2/8\pi$ \cite{Ferrer}, or Maxwell contribution. These terms are included since they also participate in the gravitational stability of the star.

Eqs.~\eqref{EoSSS} allow to analyze the influence of the magnetic field in the stability of SQM, so that
\begin{equation}\label{des}
\left.\frac{E}{n_B}\right|_{SQM}^B<\left.\frac{E}{n_B}\right|_{SQM}^{B=0}<\left.\frac{E}{n_B}\right|_{\mathrm{Fe}^{56}} = 930~\text{MeV/fm}^3,
\end{equation}
is fulfilled inside the star. Fixing $\left.E/n_B\right|_{SQM}^B=930$~MeV/fm$^3$ lead us to the stability window of SQM in the plane $\mathbf{B}$~vs~$\mathtt{B_{Bag}}$.
\begin{figure}[h]
	\centering
	\begin{tabular}{cc}
		\includegraphics[width=0.48\textwidth]{Bag_B}
		&  \includegraphics[width=0.51\textwidth]{EnbPnoB2}
	\end{tabular}
	\caption{Left panel: Stability window for magnetized SQM in the plane $\mathbf{B}$ vs $\mathtt{B_{Bag}}$. Right panel: Energy per baryon ($E/n_B$) as a function of the perpendicular pressure ($P_{\perp}$) at $B=[0, 5\times10^{17}, 10^{18}]$~G and for fixed values of the  $\mathtt{B_{Bag}}=[65, 75]$~MeV/fm$^3$.}
	\label{bagdepend}
\end{figure}
This is shown in the left panel of Fig. \ref{bagdepend}, where the SQM is stable for the pairs of $\mathbf{B}$~vs~$\mathtt{B_{Bag}}$ below the curves and unstable otherwise. We can also see from this plot that increasing the magnetic field and the baryon density augments the stability region. So, we can say that magnetic field  contributes to stabilize SQM with respect to nuclear matter in astrophysical environments, favoring the existence of SSs.
This conclusion is reinforced by the results shown in
the right panel of Fig. \ref{bagdepend}, that illustrates the energy per baryon as a function of the perpendicular pressure. At $P_{\perp}=0$, $\left.E/n_B\right|_{SQM}^B<\left.E/n_B\right|_{\mathrm{Fe}^{56}}$, i.e. the magnetic field increases the stability region, while the $\mathtt{B_{Bag}}$ diminishes it.

A more detailed stability analysis should take into account the $\mathtt{B_{Bag}}$ dependency with the magnetic field and the baryon density \cite{Ferrer}. Nevertheless, the theoretical attempts to find those dependencies or to constrain the $\mathtt{B_{Bag}}$ values are yet inconclusive and very model dependent~\cite{Buballa:2003qv}. Therefore, in what follows we take $\mathtt{B_{Bag}}$ as an independent external parameter and study SSs for two fixed, and reasonable, values of it: $65$~MeV/fm$^3$ and $75$~MeV/fm$^3$.

The SQM EoS obtained for those values of $\mathtt{B_{Bag}}$ at $B=[0, 5\times10^{17},10^{18}]$~G are depicted in Fig. \ref{eosplot}. Note that at higher values of the magnetic field, the difference between the perpendicular and parallel pressures is more appreciable. On the other hand, for a fixed $B$, the pressures decrease when increasing the $\mathtt{B_{Bag}}$.
These effects on the EoS will be reflected in the macroscopic structure of the star as we will see in the next section.
\begin{figure}[!h]
\begin{tabular}{cc}
    \includegraphics[width=0.48\textwidth]{EoSBag65}
    &\includegraphics[width=0.48\textwidth]{EoSBag75}
    \end{tabular}
    \caption{EoS for magnetized SSs at $B=[0, 5\times10^{17}, 10^{18}$~G] and for fixed values of the  $\mathtt{B_{Bag}}=[65$~MeV/fm$^3$, $75$~MeV/fm$^3$].}
    \label{eosplot}
\end{figure}
\section{Magnetized Strange Stars: Mass and Radii}\label{sec3}

The axial symmetry imposed in the star by the magnetic field is irreconcilable with the spherical symmetry of standard TOV equations. Consequently, it is desirable the use of axisymmetric metrics if one wishes to describe the structure of magnetized COs (see \cite{prd} and references therein). Here, we follow  and use a set of axisymmetric structure equations derived from the so--called $\gamma$ metric \cite{prd} \begin{align}\label{Ec(13)}
ds^2 = & - \left[1-2M(r)/r\right]^{\gamma}dt^2 + \left[1-2M(r)/r\right]^{-\gamma}dr^2
+ r^2\sin\theta d\phi^2 + r^{2}d\theta^2,
\end{align}
where $\gamma = z/r$ accounts for the deformation of the matter source with respect to the spherical shape and parametrizes the polar radius, $z$, in terms of the equatorial one, $r$.

Starting from this metric and considering the anisotropic energy-momentum tensor of magnetized matter, the following structure equations are obtained \cite{prd}
\begin{subequations}\label{gTOV}
    \begin{eqnarray}
    \frac{dM}{dr}&=& 4\pi r^{2}\frac{(E_{\parallel}+E_{\perp})}{2}\gamma, \label{gTOV1}\\
    \frac{dP_{\parallel}}{dz}&=&\frac{1}{\gamma}\frac{dP_{\parallel}}{dr}\nonumber\\
    &=&-\frac{(E_{\parallel}+P_{\parallel})[\frac{r}{2}+4\pi r^{3}P_{\parallel}-\frac{r}{2}(1-\frac{2M}{r})^{\gamma}]}{\gamma r^{2}(1-\frac{2M}{r})^{\gamma}}, \label{gTOV2}\\
    \frac{dP_{\perp}}{dr}&=&-\frac{(E_{\perp}+P_{\perp})[\frac{r}{2}+4\pi r^{3}P_{\perp}-\frac{r}{2}(1-\frac{2M}{r})^{\gamma}]}{r^{2}(1-\frac{2M}{r})^{\gamma}}. \label{gTOV3}
    \end{eqnarray}
\end{subequations}
Eqs.~(\ref{gTOV}) describe the variation of the mass and the pressures with the spatial coordinate $r$ for an anisotropic axially symmetric CO as long as the parameter $\gamma$ is close to one \cite{prd}. Note that they are coupled through the dependence with the energy density and the mass.

To solve Eqs.~(\ref{gTOV}), one starts from a point in the star center with $\mathcal E_c = E(r=0)$, $P_{\parallel_c} = P_\parallel(r=0)$ and $P_{\perp_c} = P_\perp(r=0)$, taken from the EoS, and ends the integration at $P_\parallel (Z) = 0$  and $P_\perp (R) = 0$, the conditions that defines the radii of the star. But before doing so, one should deal with the fact that $\gamma$ is a free parameter that can not be obtained from Eqs.\eqref{gTOV}. To overcome this difficulty in \cite{prd} we proposed the ansatz $\gamma = P_{\parallel_c}/P_{\perp_c}$, that connects the geometry of the system with the anisotropy produced by the magnetic field. This interpretation of $\gamma$ is based on the fact that for spherical stars a lower central pressure leads to a smaller radius for a fixed central energy, together with the relation $\gamma = z/r$. This ansatz implies that the shape of the star is only determined by the anisotropy of the EoS in its center and neglects the dependence of the deformation on the pressures inner profiles. However, its use gives reasonable results as long as $\gamma \simeq 1$  \cite{prd}. In addition, when setting $B=0$, $P_{\perp}=P_{\parallel}$, $\gamma=1$ and the standard non-magnetized solution for the structure of COs (TOV equations) is recovered.

In Fig.~\ref{gammaMR} we show the solutions of Eqs.~\eqref{gTOV} with the use of Eqs.~\eqref{EoSSS} for several values of the magnetic field and the bag energy. They are compared with the non-magnetized case and with the TOV solutions considering the pairs ($E,P_{\perp}$) and ($E,P_{\parallel}$) as independent EoS. The maximum masses and corresponding radii are shown in Table \ref{mrmaxg}. In the case of TOV solutions, using one EoS or the other leads to different mass-radius relations, whose differences increase with the magnetic field. For a given energy density range, a higher pressure implies bigger and massive stars.  Also, for a fixed value of the magnetic field and the $\mathtt{B_{Bag}}$, the difference in the stars size obtained with the pairs ($E,P_{\perp}$) and ($E,P_{\parallel}$), is larger for heavier stars. This suggest that more massive stars will have a greater deformation.

Unlike TOV equations, Eqs.~\eqref{gTOV} allow us to model the star as a spheroidal with an equatorial radius $R$ and a polar radius $Z$. So, in Fig.~\ref{gammaMR} the $M-R$ and $M-Z$ curves correspond to an unique sequence of  stars, while the $M-R_{\perp}$ and the $M-R_{\parallel}$ ones stand for two different sequences with EoS ($E,P_{\perp}$) and ($E,P_{\parallel}$) respectively.
\begin{figure}[h!]
	\centering
	\begin{tabular}{cc}
		\includegraphics[width=0.5\textwidth]{MRBag65gammaE}
		& \includegraphics[width=0.5\textwidth]{MRBag75gammaE}
	\end{tabular}
	\caption{Solutions for spheroidal configurations in comparison with TOV solutions and the non-magnetized configuration at $B=5\times10^{17}$~G and $B=10^{18}$~G, where $r$ represents both radii $R$ and $Z$. Left panel: $\mathtt{B_{bag}}=65$~MeV/fm$^3$. Right panel: $\mathtt{B_{bag}}=75$~MeV/fm$^3$~\cite{Perez2019}.}
	\label{gammaMR}
\end{figure}
The stellar configurations obtained with Eq.~\eqref{gTOV} are oblate objects ($R>Z$), as expected since $P_{\perp}>P_{\parallel}$ (see Figs. \ref{eosplot}, \ref {gammaMR}). Besides, on the contrary of what happens with TOV solutions, for which the difference between $R_{\perp}$ and $R_{\parallel}$ increases with the mass, the deformation of our spheroidal stars -the differences between the equatorial and the polar radius- decreases with the mass \cite{Perez2019} (See Fig.~\ref{gammaMR}). Hence, the importance of building a model, as the one we present, that takes into account both pressures simultaneously.
\begin{table*}[ht!]
	\caption{Maximum values of masses ($M$) Eq.~\eqref{gTOV} and the corresponding radius $R$, $Z$)}
	\label{mrmaxg}
	\centering
	\begin{tabular}{|l|l|l|l|l|}
		\hline
		$B$(G) & $\mathtt{B_{Bag}}$(MeV$/$fm$^3$) & $M(M_{\odot})$ & R(km) & Z(km)\\
		\hline
		\hline
		$5\times 10^{17}$& 65 & 1.74 & 9.79 & 9.41 \\
		& 75 & 1.63 & 9.08 & 8.79 \\

		$10^{18}$& 65 & 1.84 & 12.01 & 7.66\\
		& 75 & 1.69 & 10.20 & 7.88 \\
		\hline
	\end{tabular}
\end{table*}

The atypical behavior of the curves at $B=10^{18}$~G Fig.~\ref{gammaMR} is because, for this value of the magnetic field the $\gamma$ parameter escapes the range established in \cite{prd}, ($1>\gamma>0.8$) for the validity of Eqs.~\eqref{gTOV}. It is interesting that other models of Neutron Stars EoS and structure equations also present difficulties when treating magnetic fields of the order or higher than $B=10^{18} $~G. For example, in \cite{Daryel2015}, where the anisotropy is considered through a metric in cylindrical coordinates, the metric coefficients diverge for $B \geq 1.8\times 10^{18} $~G. Besides, the theoretical limit established by the Virial theorem for the maximum magnetic field that can be supported by a NS is precisely $B\simeq10^{18}$~G~\cite{Khar}. Therefore, our results seem to support this limit, indicating the relevance of going deeper in the understanding of the physical reasons behind it.
\begin{figure}[!h]
	\includegraphics[width=0.48\textwidth]{g_E_gamma}
	\caption{$\gamma$ parameter as a function of central energy density at $B=[0, 5\times10^{17}, 10^{18}$~G] and for fixed values of the  $\mathtt{B_{Bag}}=[65$~MeV/fm$^3$, $75$~MeV/fm$^3$].}    \label{gammaE}
\end{figure}

The effects of varying $\mathtt{B_{Bag}}$ and the magnetic field on the deformation of the stars can also be seen through the ellipticity, defined as \cite{Rizaldy}
\begin{equation}\label{elip}
\epsilon=\sqrt{1-\gamma}.
\end{equation}
In the spherical case  $\gamma\rightarrow 1$ and $\epsilon\rightarrow 0$, while for the most deformed stars, $\gamma\rightarrow 0$ and $\epsilon\rightarrow 1$. In Fig.~\ref{emgamma} we show the ellipticity $\epsilon$ as a function of the  mass and the central energy density, and show that the increase of the magnetic field and the $\mathtt{B_{Bag}}$ leads to a greater deformation of the star.
\begin{figure}[h!]
	\centering
	\begin{tabular}{cc}
		\includegraphics[width=0.5\textwidth]{Eccent_M} &
	    \includegraphics[width=0.5\textwidth]{Eccent_M_E}
	\end{tabular}
	\caption{Elipticity ($\epsilon$) as a function of the mass $M$ (left) and the energy density  $E$ (right) at $B=[0, 5\times10^{17}, 10^{18}$~G] and for fixed values of the  $\mathtt{B_{Bag}}=[65$~MeV/fm$^3$, $75$~MeV/fm$^3$].}
	\label{emgamma}
\end{figure}
\subsection{Stability of the solutions of the $\gamma$--equations}

The mass-radii curves obtained for SSs with Eqs.~(\ref{gTOV}) correspond to stellar configurations in hydrostatic equilibrium \cite{Glendenning}. However, equilibrium does not necessarily imply stability~\cite{Shapiro, Gleiser1988}. To find which of those solutions represents stable COs we will use two criteria. The first is related to the stability of the star with respect to radial oscillations and requires that $dM/dE_0>0$ \cite {Shapiro}, so the CO will not fall apart. The second requires the star's gravitational mass (the one calculated with the structure equations) to be lower than its baryonic mass ($M_B$), which is  the sum of the masses of all its particles. This last criterion is applied to guarantee the stability of the CO with respect to the dispersion of the matter that composes it~\cite{Gleiser1988}.

To analyze the stability with respect to radial oscillations of the $\gamma$-equations solutions, in the left panel of Fig.~\ref{megamma} the mass has been plotted as a function of the central energy density. The points represent the maximum of each curve. Therefore, the stars to the left of the maximum are stable ($dM/dE_0>0$), while the ones to the right are not ($dM/dE_0<0$) and must be discarded. For fixed magnetic field, increasing $\mathtt{B_{Bag}}$ increase the stability region, whereas, for a fixed $\mathtt{B_ {bag}}$, increasing $B$ reduces it. Note that the solutions for $B=10^{18}$~G are out of the validity range of $\gamma$--equations for low densities while are unstable with respect to radial oscillations for the higher ones. Therefore, in what follows, we will not consider the observables for this value of the magnetic field.
\begin{figure}[h!]
	\centering
	\begin{tabular}{cc}
		\includegraphics[width=0.5\textwidth]{M_E_gammaE}
		& \includegraphics[width=0.5\textwidth]{MMb_E_gamma}
	\end{tabular}
	\caption{Mass $M$ as function of energy density  $E$ (left) and ratio between $M$ and the baryon mass $M_B$ (right) at $B=[0, 5\times10^{17}, 10^{18}$~G] and for fixed values of the  $\mathtt{B_{Bag}}=[65$~MeV/fm$^3$, $75$~MeV/fm$^3$].}
	\label{megamma}
\end{figure}

To apply the second stability criterion we start from the definition of the baryonic mass $M_B$\cite{Mauro}
\begin{equation}
M_B=m_N\int_{0}^{R}\frac{4\pi r^2 n_B(r)}{[1-2Gm(r)/r]^{1/2}}dr,
\end{equation}
where $m_N$ is the neutron mass and $n_B(r)$ is the baryonic density. The regions where $ M/M_B<1 $ are  stable, while those where $M/M_B>1$ correspond  to unstable stellar configurations. Therefore, these regions should also be discarded. But, as shown in the right panel of Fig. \ref{megamma}, all the stellar configurations here obtained are stable with respect to this stability criterion.

\subsection{Comparison of the model with Strange Stars candidates}

In this subsection, we compare the spheroidal static configurations obtained from $\gamma$-equations with some observable data for candidates to be SSs (See Table \ref{candidatas}). As spheroidal SSs have two main radii, for the comparison we defined a mean radius $R_m$, so that the sphere it determines is equal to the surface of the spheroidal star
\begin{equation}
A=2\pi R\left[R+\frac{Z}{\epsilon}\arcsin\epsilon\right],
\end{equation}
where $\epsilon$ is the ellipticity. In this way, the radius $R_m$ could be connected with the surface of the COs and consequently, with its emission properties\cite{Debarati}.

Fig.\ref{mrcang} shows our theoretical results for the masses and radii together with the observational ones for the candidates. The shaded regions correspond to regions forbidden by theoretical requirements (see the legend of Fig.\ref{mrcang}). The curves at $B=0$ are in the allowed region. In the case of the curves at $B=5\times10^{17}$~G, only the one corresponding to $\mathtt{B_{Bag}}=65$~MeV$/$fm$^3$ has a small  section in the region of rotational instability. In general, our curves are in the region where the candidates are concentrated. In particular, the values of the mass and radius of the stars Vela X-1 and 4U1608-52 fall on the curve $B = 5\times 10^{17}$~G for $\mathtt{B_ {bag}}=65$~MeV$/$fm$^3$. This supports our model of magnetized Strange Stars, meaning it is consistent with the observations.
\begin{figure}[h!]
	\centering
	\begin{tabular}{cc}
		\includegraphics[width=0.5\textwidth]{MRmCan65.pdf}
		& \includegraphics[width=0.5\textwidth]{MRmCan75.pdf}
	\end{tabular}
	\caption{Comparison of observational data for candidates of SSs with the stars obtained from Eqs.~\eqref{EoSSS} and \eqref{gTOV}. The shaded regions correspond to theoretical constraints. Gravitational collapse (I). Requirement of finite pressure inside the star (II). Causality (III). Rotational stability (IV) \cite{Lattimer_prog}. Left panel: $\mathtt{B_{Bag}}=65$~MeV/fm$^3$. Right panel: $\mathtt{B_{Bag}}=75$~MeV/fm$^3$.}
\label{mrcang}
\end{figure}
\section{Magnetized Strange Stars: Redshift and Mass Quadrupole}\label{sec4}

\subsection{Gravitational Redshift}

Gravitational redshift is one of the main effects predicted by the Theory of General Relativity and in turn, constitutes one of its fundamental tests. The frequency of an atomic clock depends on the value of the gravitational potential in the place of its location. So when a photon is observed from a point at a higher gravitational potential, its wavelength is redshifted.

For a spheroidal static CO, the redshift is given by \cite{Zubairi2015}
\begin{equation}
z_{rs}=\frac{1}{\big(1-\frac{2M}{R}\big)^{\gamma/2}}-1,\label{zgamma}
\end{equation}
where $M$ and $R$ are, respectively, the mass and the equatorial radius of the CO. If $\gamma = 1$ the spherical case is recovered.

The definition of $z_{rs}$ explicitly includes the term $M/R$, so its measurement can narrow these parameters. In addition, $z_{rs}$ has a maximum at the star's maximum point of mass, so it can also be used to rule out EoS that do not lead to observable redshifts and therefore constitutes a benchmark to evaluate the feasibility of NSs models \cite{Camenzind}. Our results for $z_{rs}$ are shown in Fig.~\ref{zmtov} for TOV and  $\gamma$-equations solutions.
\begin{figure}[h]
	\centering
	\begin{tabular}{cc}
		\includegraphics[width=0.5\textwidth]{ZMall65}
		& \includegraphics[width=0.5\textwidth]{ZMall75}
	\end{tabular}
	\caption{Gravitacional redshift ($z_{rs}$) as a function of the mass ($M$) at $B=[0, 5\times10^{17}, 10^{18}$~G], for TOV and Eqs. $\gamma$. Left panel:  $\mathtt{B_{Bag}}=65$~MeV$/$fm$^3$. Right panel:  $\mathtt{B_{Bag}}=75$~MeV$/$fm$^3$.}
	\label{zmtov}
\end{figure}

In the spherical case, the higher values for $z_{rs}$ are obtained for the pairs ($E, P_{\perp}$). In this case, for a fixed $\mathtt{B_{Bag}}$, increasing the magnetic field increases $z_{rs}$, while changing $\mathtt{B_{Bag}}$ barely affects it. If we compare the results for $B=0$ with the solutions of $\gamma$-equations for $B = 5\times 10^{17}$~G, we see that with these, lower values of $z_{rs}$ are obtained and the curves are far apart. This difference can be very useful to discriminate between models with/without magnetic field  when comparing with observational values.

\subsection{Mass quadrupole moment}

Finally, we calculate the mass quadrupole moment of the spheroidal stars. This magnitude is directly related to the amplitude of the GWs \cite{Magdalena}, since they are only emitted in situations where an asymmetry of mass is generated that gives rise to a quadrupolar moment. Therefore, spherical stars do not have a quadrupole moment and cannot generate GWs. In contrast, magnetized stars, being deformed, have non-zero quadrupole moments that will contribute to their GWs emission. In the framework of our structure equations, the quadrupole moment of the SSs is \cite {Herrera99}
\begin{equation}\label{Mq}
Q=\frac{\gamma^3}{3}M^3(1-\gamma^2).
\end{equation}
where for $\gamma =1$, $Q=0$, as corresponds to the spherical case.

\begin{figure}[h!]
	\centering
	\begin{tabular}{c}
		\includegraphics[width=0.5\textwidth]{M_Q_gamma}
	\end{tabular}
	\caption{Mass quadrupole moment ($Q$) as a function of the mass ($M$) at $B=5\times10^{17}$~G for $\mathtt{B_{Bag}}=65$~MeV$/$fm$^3$ and $\mathtt{B_{Bag}}=75$~MeV$/$fm$^3$.}
	\label{qmgamma}
\end{figure}

Fig.\ref{qmgamma} shows the SSs quadrupole moment as a function of the star mass. The oscillations in the curve are an effect of the presence of the sum by the Landau levels in the EoS. $Q$ diminishes with $\mathtt{B_{Bag}}$ and its maximum is reached for stars in the region of intermediate masses and deformation. This behavior is due to the simultaneous dependence of $Q$ on $M$ and $\gamma$ and, in particular, is determined by the fact that $\gamma$ depends on the EoS and therefore varies between the stars. This result is different from the one obtained in \cite{Zubairi2015}, where structure equations derived from the $\gamma$ metric were solved by taking $\gamma$ as a free parameter set by hand. In that case, since the mass quadrupole is leaded by the mass, its highest values are reached for the more massive stars. Therefore, connecting $\gamma$ with the physics of the problem has a direct impact on the observables, that again could serve as a way to discriminate between models.

\section{Conclusions}\label{sec5}

This work is an step forward in our studies about magnetized SSs \cite{Felipe2008,Felipe2009JPhG, ARD}. We re-analyzed the stability of SQM under the action of a magnetic field and found that its presence reinforces the Bodmer-Witten's conjecture. To compute the star's observables the EoS have been restricted to the stability region (energies below $930$~MeV). The macroscopic properties of the SSs --mass, radius, deformation, gravitational redshift, mass quadrupole moment-- were calculated using the $\gamma$-equations~\cite{prd} for spheroidal compact objects. The mass--radius curves of the stable configurations obtained are consistent with the observed properties of SSs candidates and comply with the theoretical constraints; this result support the plausibility of our model.

In our model, less massive stars suffer bigger deformations in contrast with the results from TOV solutions for the perpendicular and the parallel pressure independently. This reveals the model--dependency of the results and highlights how important is the construction of even more realistic models. On the other hand, augmenting both, $\mathtt{B_{Bag}}$ or $B$, increases the deformation. However, the effects of changing the bag energy and the magnetic field in the stability with respect to radial oscillations of the solutions of $\gamma$--equations opposes: the stability region increases with $\mathtt{B_{Bag}}$ and decreases with $B$. Nevertheless, none of these parameters influences the stability of the star with respect to the baryon mass criterion.

Another interesting feature of our model is that the redshift and the mass quadrupole moment depend explicitly on the deformation through the EoS --because $\gamma$ appears on their mathematical expressions. As a consequence, the curve of the  gravitational redshift of SSs has remarkable differences with respect to the ones of the spherical case. Besides, the maximum values of the mass quadrupole moment occur for the stars with intermediate values of masses and deformation, suggesting that these are the stars that should produce the most intense GWs emission, instead of those that have the maximum masses or deformation. This result is as interesting as unexpected and deserves a deeper research in a future work.

\section*{Acknowledgments}

The authors have been supported by the grant No. 500.03401 of the PNCB-MES Cuba and by the ICTP Office of External Activities through NT-09.